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SUMMARY 


This  study  reports  on  observations  of  dawn/dusk  transition  fading  of  very 
low  frequency  radio  waves  on  a  long,  transequatorial  northeasterly  propagation 
path.  A  mode  conversion  model  is  developed  and  numerical  modeling  results  are 
presented  for  transition  fading  on  long  transequatorial  paths,  both  northeast¬ 
erly  and  southwesterly.  From  these  results  the  following  general  conclusions 
are  drawn: 

a.  Anomalous  transequatorial  transition  fading  is  not  observed  for  the 
northeasterly  propagation  path.  The  transition  fading  data  support  completely 
the  mode  conversion  model  of  transition  fading  developed  by  Crombie. 

b.  The  transequatorial  sunrise  transition  reported  in  the  literature  as 
anomalous  is  completely  explained  by  correctly  allowing  for  variations  of  the 
geomagnetic  field  along  the  path.  The  geomagnetic  field  influences  the  night¬ 
time  modal  parameters  to  such  an  extent  for  the  southwesterly  path  considered 
that  the  approximate  WKM  method  of  mode  summation  cannot  be  employed. 

c.  Careful  application  of  modal  techniques  appears  to  be  capable  of 
describing  all  LF  and  VLF  propagation  in  the  earth-ionosphere  waveguide, 
provided  the  effect  of  the  earth's  magnetic  field  is  properly  included  at 
night.  However,  an  improved  model  of  the  ionosphere  at  altitudes  below  100  km 
is  desired. 
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INTRODUCTION 
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The  propagation  of  low  and  very  low  frequency  ( LF  and  VLF )  radio  waves 
has  been  of  interest  since  the  beginning  days  of  radio.  The  epoch-marking 
experiments  at  the  beginning  of  this  century  by  Marconi  were  made  using  what 
is  now  called  the  LF  band.  These  experiments  gave  impetus  to  the  development 
of  LF  and  VLF  radio  communications  and  to  studies  of  their  propagation  charac¬ 
teristics,  particularly  to  the  variation  of  signal  strength  with  distance  from 
the  transmitter. 

There  are  many  practical  disadvantages  to  using  these  frequency  bands 
(Watt,  1967).  These  include  the  very  high  atmospheric  noise  level,  which 
necessitates  high  power  transmissions,  the  limited  bandwidth  available  and  the 
high  cost  and  low  efficiency  of  transmitting  antennas.  These  disadvantages 
did  not  exist  at  shorter  wavelengths  and  the  development  of  shortwave  radio  in 
the  1920s  led  to  the  near  demise  of  interest  in  the  longer  wavelengths. 

For  a  variety  of  reasons,  interest  in  VLF  radiowave  propagation  was 
revived  and  has  greatly  increased  during  the  last  two  decades.  Near  the 
earth's  surface,  the  structure  of  the  VLF  electromagnetic  field  depends  on  the 
properties  not  only  of  the  earth  but  also  of  the  lower  regions  of  the  ionos¬ 
phere  between  50  to  100  km,  the  D-region  and  lower  E-region.  VLF  signals  are 
not  greatly  affected  by  most  ionospheric  disturbances,  apart  from  polar  cap 
events  (Al'pert  and  Fligel,  1970),  and  the  use  of  this  frequency  band  affords 
a  very  useful  tool  for  studying  the  upper  regions  of  the  ionosphere,  out  to 
several  earth  radii.  This  use  stems  from  the  fact  that  VLF  signals,  either 
man-made  or  natural,  from  electric  storms  fo>  example,  propagate  through  the 
lower  regions  of  the  ionosphere  and  are  guided  by  ducts  aligned  along  the  geo¬ 
magnetic  field  lines  (the  so-called  whistler  mode)  whereupon  they  enter  the 
earth's  other  hemisphere.  Thus  analysis  of  such  signals  can  yield  information 
about  the  near  earth  plasma  and  in  fact  provided  the  first  evidence  for  the 
existence  of  a  plasmasphere  above  the  F-region.  Furthermore,  VLF  has  been 
found  to  be  perhaps  the  most  sensitive  terrestrial  means  of  detecting 
ionospheric  disturbances  produced  by  the  entry  of  solar  particles  into  the 
ionosphere  (potemra  and  Rosenberg,  1973). 
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Practical  application  of  VLF,  particularly  from  10  to  30  kHz  has  also 
stimulated 'tinterest  in  this  frequency  band  -  VLF  signals  exhibit  low  attenua¬ 
tion  (on  the  order  of  a  few  dB  per  1000  km) ,  and  this  frequency  band  is  widely 
used  in  low  data  rate,  global  communication  systems.  The  phase  of  VLF  trans¬ 
mission  is  highly  stable  and  undergoes  nearly  reproducible  diurnal  variations. 
It  has  been  found,  for  example  (Belrose,  1968),  that  the  stability  of  VLF 
transmissions  is  sufficient  to  permit  frequency  comparisons  to  within  a  few 
parts  in  10  ^  which  is  a  few  powers  of  ten  better  than  is  possible  at  high 
frequencies.  This  phase  stability  leads  to  applications  in  navigational  sys¬ 
tems,  frequency  comparisons  and  clock  synchronization  or  timing  systems. 

The  ability  to  determine,  through  VLF  propagation  theory,  the  variations 
of  VLF  radiowaves  as  a  function  of  distance  is  of  importance  to  the  navigation 
or  communications  engineer,  who  must  design  his  system  for  optimum  performance 
in  his  total  operational  environment,  which  includes  the  propagation  environ¬ 
ment.  Such  theoretical  calculations  require  an  adequate  description  of  the 
propagation  environment.  A  comparison  of  calculated  propagation  parameters 
with  those  measured  experimentally  enables  information  to  be  obtained 
regarding  the  parameters  of  the  propagation  medium.  A  knowledge  of  these 
parameters  is  important  for  the  study  of  physical  processes  which  occur  in  the 
propagation  medium. 

A.  VLF  Radiowave  Propagation  to  Great  Distances 

1 .  Theoretical  Developments 

This  study  is  concerned  with  some  aspects  of  the  propagation  of  VLF 
radiowaves  to  great  distances,  i.e.,  with  their  amplitude,  polarization  and 
phase  variations.  VLF  propagation  to  great  distances  has  been  considered  in  a 
large  number  of  theoretical  papers,  using  a  variety  of  formulations.  We  will 
not  attempt  to  review  this  massive  literature  but  will  concentrate  on  some 
general  ideas  that  have  been  incorporated  into  the  various  formulations. 

One  of  the  major  formulations  for  consideration  of  VLF  propagation  to 
great  distances,  the  VLF  wave-hop  theory,  has  been  presented  by  Berry  and 
Chrisman  (1965) .  In  this  formulation,  the  solution  for  propagation  of  VLF 
waves  between  a  spherical  earth  and  a  concentric  ionosphere  is  developed  in  a 
rigorous  complex  integral  representation  and  the  integrand  is  then  expanded  in 
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terms  of  a  geometrical-type  series.  With  the  order  of  summation  and  integra¬ 
tion  exchanged,  evaluation  of  the  complex  integrals  using  saddle  point  approx¬ 
imations  leads  to  the  identification  of  the  series  as  the  ray-hop  series  of 
geometric  optics.  The  complex  integrals  are  thus  called  wave-hops.  For  dis¬ 
tances  near  or  beyond  the  caustic,  numerical  integration  or  application  of 
residue  theory  must  be  employed  to  evalute  the  integrals.  The  wave-hop  formu¬ 
lation  requires  both  azimuthal  and  longitudinal  invariance  of  the  earth  and 
ionosphere,  a  situation  which  exists  only  for  short  ranges. 

Johler  (1970)  has  presented  a  theory  for  radiowave  propagation  from  VLF 
to  MF  in  terms  of  spherical  wave  functions.  As  with  the  wave-hop  formulation, 
the  spherical  wave  formulation  is  not  applicable  to  consideration  of  varia¬ 
tions  of  the  earth  or  ionosphere  in  either  the  longitudinal  or  azimuthal 
directions.  The  spherical  wave  function  formulation  has  not  been  widely  used 
because  of  its  complexity  and  difficulty  of  numerical  implementation. 

VLF  propagation  to  great  distances  can  be  conveniently  represented  in 
terms  of  waveguide  mode  propagation,  where  the  finitely  conducting  curved 
earth  and  the  anisotropic,  imperfectly  conducting  curved  ionosphere  with 
dipping  magnetic  field  form  the  boundaries  of  a  waveguide.  The  basic  idea 
that  the  earth  and  lower  ionosphere  form  a  waveguidi  actually  goes  back  to  the 
work  of  Watsc  (1918,  1919),  who  postulated  that  the  Kennelly-Heaviside  layer 

c°uld  be  represen' ed  i.s  a  conducting  shell  concentric  with  the  spherical  con¬ 
ducting  earth.  Watson  <1  ,v«  loped  an  exact  slowly  convergent,  harmonic  series- 
type  solution  to  the  F  propagation  problem  and  further  devised  a  transforma¬ 
tion  of  the  series  to  a  highly  convergent  form.  The  transformation,  which 
involves  distortion  of  the  contour  of  integrat  ion  in  the  complex  wave  number 
plane,  is  now  identified  with  Watson's  name,  although  the  transformed  series 
is  often  ’rediscovered"  by  subsequent  workers.  Although  "wavegu i des"  were  not 
really  "discovered"  until  the  1930s,  tire  terms  in  Watson's  ti ansf ormed  series 
were  simply  waveguide  modes. 

The  theory  of  VLF  waveguide  mode  propagation  has  been  presented  in  the 
literature  in  two  different  formulation-,  one  formulation,  based  on  the 
pioneering  work  of  Watson  (19  19)  and  more  recently  in  the  theory  given  by 
Schumann  (1952,  1954)  and  developed  extensively  by  Wait  (1962),  treats  the 
problem  in  spherical  coordinates  wherein  the  fields  are  expanded  in  terms  of 


azimuthal,  longitudinal  and  radial  functions.  In  a  sequence  of  simplifying 
approximations,  including  the  assumption  of  azimuthal  and  longitudinal  inde¬ 
pendence  of  both  ground  and  ionospheric  properties.  Wait  (1962)  develops  a 
modal  (eigenvalue)  equation  in  terms  of  Airy  functions.  The  Airy  functions 
come  about  as  a  result  of  a  third-order  approximation  to  spherical  wave  func¬ 
tions  (Watson,  1944).  The  solution  of  the  eigenvalue  equation  yields  the  pro¬ 
pagation  characteristics  of  the  VLF  modes  in  the  earth-ionosphere  waveguide 
and  consideration  of  orthogonality  properties  of  the  radial  (Airy)  functions 
yields  excitation  factors  of  the  modes.  The  excitation  factor  is  a  quantity 
giving  the  amplitude  of  the  wave  excited  in  a  given  mode  by  a  given  source. 
Another  simplifying  approximation  often  made  involves  the  replacement  of  the 
inhomogeneous  ionosphere  with  an  impedance  boundary,  placed  at  a  height  where 
the  bulk  of  the  VLF  energy  is  assumed  to  be  reflected.  A  similar  impedance 
boundary  replacement  is  also  made  for  the  earth.  The  impedance  boundary 
formulation  permits  the  specification  of  the  ratio  of  an  electric  field 
component  to  a  magnetic  field  component  at  the  boundary  without  further  con¬ 
sideration  of  the  region  beyond  the  boundary.  Analytically,  formulation  in 
terms  of  impedance  boundaries  is  equivalent  to  assumption  of  "homogeneous 
boundary  conditions"  (Friedman,  1959;  Morse  and  Feshbach,  1953).  The  replace¬ 
ment  of  an  inhomogeneous,  isotropic  ionosphere  by  an  "equivalent"  impedance 
boundary  is  a  straightforward  procedure.  Such  a  replacement  for  an  anisotro¬ 
pic  ionosphere  is  greatly  complicated  by  the  coupling  of  transverse  magnetic 
and  transverse  electric  polarizations.  Boundary  conditions  must  then  be  form¬ 
ulated  in  terms  of  an  impedance  matrix.  Unfortunately,  the  complication  is 
sufficient  to  render  many  of  Wait's  results  inapplicable  to  the  highly  aniso¬ 
tropic  nighttime  ionosphere  and  to  limit  their  application  to  the  daytime 
ionosphere,  for  which  the  effect  of  anisotropy  is  slight.  The  theory  has  been 
significantly  modified  and  extended  to  anisotropic  ionospheres  by  Galejs 
( 1972)  . 

Another  formulation  of  the  earth-ionosphere  waveguide  problem  has  been 
presented  by  Budden  (1961a).  It  is  Budden's  formulation  that  forms  the  basis 
for  this  study  and  will  be  discussed  in  more  detail  in  following  sections.  In 
this  formulation,  a  modal  equation  is  developed  in  terms  of  reflection  coeffi¬ 
cients  of  the  ionosphere  and  the  earth.  The  formulation  in  terms  of  reflec- 
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tion  coefficients  permits  determination  of  propagation  characteristics  of 
waveguide  modes  for  essentially  any  ionospheric  or  ground  conditions,  limited 
only  by  the  ability  to  compute  reflection  coefficients  for  the  environment 
considered.  Budden's  formulation  is  essentially  planar,  although  he  does 
point  out  how  to  include  earth  curvature  in  the  direction  of  propagation  in  an 
approximate  way  by  modifying  the  refractive  index  in  the  space  between  the 
earth  and  the  ionosphere  (Budden,  1962) .  A  more  rigorous  technique  was 
introduced  by  Richter  (1966)  and  has  been  employed  in  VLF  propagation  studies 
by  Pappert  (1968)  and  Abbas  et  al.  (1975).  This  technique,  which  is  employed 
in  this  study  and  will  be  discussed  more  fully  in  following  sections,  is  to 
assume  cylindrical  stratification  and  effect  a  conformal  transformation  from 
cylindrical  to  Cartesian  coordinates.  Modal  excitation  values  are  determined 
in  Budden's  formulation  by  application  of  residue  theory  rather  than  using  the 
modal  orthogonality  properties  of  conventional  waveguide  theory  employed  by 
Wait. 


Whichever  formulation  is  used  at  the  outset,  the  use  of  waveguide  mode 
theory  leads  to  specification  of  the  VLF  field  as  a  summation  of  terms  called 
modes.  Because  most  VLF  transmitters  in  common  use  radiate  a  vertically 
polarized  field,  only  the  radial  or  vertical  component  of  the  electric  field 
is  usually  considered.  For  a  time  harmonic  source,  the  VLF  mode  sum  for  the 
vertical  electric  field  may  be  written  as 


E  (d) 


K(P,f)  ,  T  R  _lkoSnd 

— - - —  A  G  G  e 

/ — ; - TT  n  n  n 

» sin  d/a  n 


(1) 


where  K(P,f)  is  a  complex  constant  dependent  on  transmitted  power  (p)  and 
frequency  (f);  d  is  the  distance  from  the  transmitter  on  a  homogeneous,  smooth 
earth  of  radius  a;  is  the  excitation  factor  for  mode  n,  normaiized  to  unity 
for  flat  earth,  perfectly  conducting  boundaries;  kQ  is  the  free  space  wave 
number;  Sn  is  the  propagation  factor;  and  G^,R  represents  height-gain  func¬ 
tions  for  mode  n,  normalized  to  unity  at  the  ground.  One  height-gain  function 
is  needed  for  the  transmitter  (T)  and  one  for  the  receiver  (R).  Generally, 

both  An  and  Sn  are  complex.  The  real  part  of  Sn  determines  the  distance 
dependence  of  the  phase  for  a  mode  while  the  imaginary  part  of  Sn  determines 
the  attenuation  rate.  Thus 
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(2) 


a  =  -k  Im(S  ) 
n  on 

v  /C  =  1/Re( S  ) 
n  n 

where  an  is  the  modal  attenuation  rate,  vn  is  the  phase  velocity  and  c  is  the 
speed  of  light  in  free  space. 

Numerous  assumptions  and  approximations,  some  of  which  will  be  removed  in 
following  sections,  are  implicit  in  equation  (1).  It  is  assumed  that  the 
propagation  environment  is  homogeneous  in  all  but  the  vertical  direction. 
Thus,  the  propagation  factor  Sn  is  independent  of  position.  The  influence  of 
the  curvature  of  the  earth  transverse  to  the  propagation  direction  is  approxi¬ 
mately  accounted  for  in  the  (sin  d/a)  term.  The  energy  flux  for  a  wave  tra¬ 

veling  outwards  with  cylindrical  symmetry  on  a  flat  earth  without  attenuation 
is  proportional  to  1/r.  For  a  spherically  curved  earth  this  energy  flux  would 
be  proportional  to  1/sin(r/a).  Thus  the  usual  r- term  for  a  cylindrical 

wave  is  replaced  by  the  ( sin ( r/a ) ) - ^ ^  term.  A  careful  derivation  of  the  mode 

equation  in  spherical  coordinates,  such  as  that  presented  by  Wait  (1962)  and 
discussed  previously,  yields  the  sin(d/a)  term  as  a  consequence  of  an  asymp¬ 
totic  approximation  to  a  Legendre  function.  The  mode  sum  in  equation  (1) 
further  assumes  that  "round-the-world"  signals  are  absent,  an  approximation 
not  valid  when  receiver  locations  are  near  the  antipode  of  the  transmitter; 
multi-path  focusing  then  occurs. 

Extensive  results  of  numerical  calculations  of  VLF  waveguide  mode  con¬ 
stants  have  been  presented  in  the  literature.  Early  results  given  by  Wait  and 
Spies  (1964)  show  some  general  properties  of  VLF  modes  for  typical  daytime  and 
nighttime  ionospheres  with  various  ground  conductivities  and  for  propagation 
in  the  east  to  west  or  west  to  east  directions  at  the  magnetic  equator  or  when 
the  earth's  magnetic  field  is  ignored.  in  general,  the  attenuation  rate  of 
the  lowest  order  mode  is  less  at  night  than  during  day,  and  is  less  for  west 
to  east  propagation  than  for  east  to  west.  The  attenuation  rate  increases  as 
ground  conductivity  decreases  and  does  so  rapidly  at  low  conductivities  such 
as  occur  for  the  polar  icecaps.  Each  successive  mode  order  suffers  a  larger 
attenuation  rate  than  the  previous  order  mode  for  both  day  and  night. 

Typically,  daytime  attenuation  rates  for  the  first  order  mode  range  from 
2  to  4  dB/1000  km  over  the  VLF  band  for  sea  water  ground  conditions  (essen- 
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tially  equivalent  to  a  perfect  conductor  at  VLF).  The  nighttime  attenuation 
rates  range  from  about  0.5  to  2  dB/1000  tan  over  the  VLF  band  for  the  same 
ground  conductivities. 

The  results  of  Waif  and  Spies  show  that  the  modal  phase  velocities  for 
the  curved  earth  can  be  less  than  the  free  space  speed  of  light  at  intermedi¬ 
ate  to  higher  VLF  frequencies  but  increase  markedly  at  lower  frequencies. 
Moreover,  phase  velocities  increase  at  low  ionosphere  heights  (daytime)  com¬ 
pared  to  higher  heights  (night)  and  a  decrease  in  ground  conductivity  produces 
a  decrease  in  phase  velocity.  An  increase  in  phase  velocity  with  mode  number 
is  also  indicated. 

The  modal  excitation  factors  and  height  gains  given  by  Wait  and  Spies 
show  little  dependence  on  the  terrestrial  magnetic  field.  At  higher  frequen¬ 
cies  for  isotropic  conditions,  the  excitation  factor  for  the  first  mode  is 
significantly  reduced  in  comparison  with  that  for  a  flat  earth  case.  The 
second  mode  exhibits  a  much  smaller  reduction  at  higher  frequencies.  Gen¬ 
erally,  lower  ground  conductivities  produce  an  increase  in  excitation  fac¬ 
tor.  The  excitation  factor  for  the  second  order  mode  is  significantly  higher 
than  that  for  the  first  order  mode,  especially  at  higher  frequencies  and  at 
night.  At  night  the  second  order  mode  may  be  dominant  at  the  ground  to  dis¬ 
tances  of  several  thousand  km.  For  the  first  mode  the  effect  of  height-gain 
is  to  cause  the  dominance  of  the  second  mode  to  be  greatly  diminished  or  even 
reversed.  This  height-gain  increase  with  height  for  the  first  order  mode  has 
been  called  the  whispering  gallery  effect  by  Budden  and  Martin  (1962)  or  the 
earth-detached  mode  by  Wait  and  Spies  (1963).  It  involves  successive  reflec¬ 
tions  from  the  ionosphere  with  essentially  no  intermediate  bounce  at  the 
ground. 

The  modal  constants  characteristics  presented  above  become  considerably 
more  complicated  with  full  consideration  of  the  geomagnetic  field.  From 
earliest  days  of  VLF,  experimental  investigation  provided  evidence  of  an 
apparent  violation  of  the  reciprocity  principle.  Found  (1925),  for  example, 
found  that  VLF  radiowaves  apparently  suffered  significantly  larger  attenuation 
for  propagation  in  a  generally  easterly  direction  than  in  the  reverse  direc¬ 
tion.  Some  very  early  results  presented  by  Wait  (1961)  for  sharply  bounded 
model  ionospheres  suggested  that  simple  harmonic  functions  may  adequately 
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describe  the  azimuthal  dependence  of  mode  constants  under  daytime  conditions. 
Further,  a  numerical  modeling  study  by  Ferguson  (1968)  showed  that,  for  day¬ 
time  conditions,  only  the  horizontal  component  of  the  geomagnetic  field  which 
is  transverse  to  the  direction  of  propagation  is  important.  This  component  of 
the  magnetic  field  is  the  source  of  nonreciprocity,  and  the  remaining  compo¬ 
nents  were  commonly  assumed  only  to  alter  the  propagation  characteristics 
somewhat  (Makorov  et  al.  ,  1970).  The  most  severe  geomagnetic  field  effects 
were  thus  expected  in  the  equatorial  region  for  propagation  across  the  mag¬ 
netic  meridians. 

Perhaps  the  first  indications  of  extreme  complexity  in  modal  behavior  due 
to  the  geomagnetic  field  are  found  in  the  numerical  modeling  results  of  Snyder 
( 1968a, b)  and  Pappert  (1968).  Snyder  (1968a)  addressed  the  problem  of  mode 
numbering  for  typical  nighttime  ionospheres  as  a  function  of  azimuthal  varia¬ 
tion  for  a  horizontal  magnetic  field.  Results  indicated  that  a  mode  with 
lowest  phase  velocity  in  the  west  to  east  direction  evolved  as  a  result  of 
continuous  variation  in  magnetic  azimuth  into  a  mode  in  the  reverse  direction 
with  only  second  lowest  phase  velocity.  Thus,  numbering  modes  in  terms  of 
increasing  phase  velocity,  as  is  commonly  done,  results  in  the  first  order 
mode  in  the  west  to  east  direction  evolving  into  the  second  order  mode  in  the 
opposite  direction  and  vice  versa.  Further  analysis  (Snyder,  1968b)  indicated 
that,  for  a  certain  inclination  of  a  dipping  geomagnetic  field,  and  at  a 
specific  azimuth  (near  the  north-south  direction),  the  two  modes  considered 
above  became  degenerate,  forming  a  single  mode.  For  inclinations  more  nearly 
horizontal  than  the  degenerate  conditions,  mode  numbering  inconsistencies 
resulted  but,  for  more  nearly  vertical  inclinations,  inconsistencies  in  mode 
numbering  vanished.  Pappert  (1968)  presented  the  results  of  numerical  model¬ 
ing  for  easterly  propagation  at  midlatitudes.  Assuming  a  vertical  dipole 
exciter,  Pappert  found  that,  at  the  high  frequency  end  of  the  VLF  band  (only 
near  30  kHz),  modes  which  are  principally  transverse  electric  (TE)  may  be  of 
importance  in  a  mode  sum.  Note  that  a  TE  mode  can  be  excited  by  a  vertical 
dipole  only  because  of  the  presence  of  the  geomagnetic  field. 

Snyder  and  Pappert  (1969)  extended  the  analysis  to  consider  both  easterly 
and  westerly  midlatitude  nighttime  propagation  throughout  the  VLF  band  as  well 
as  to  include  azimuthal  dependencies  of  mode  parameters  for  a  central  VLF  fre- 
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quency  at  both  midlatitudes  and  equatorial  latitudes.  It  was  found  that  the 
importance  of  principally  TE  modes  is  much  more  pronounced  for  westerly  propa¬ 
gation  at  midlatitudes  than  for  easterly,  and  that  the  influence  of  the  prin¬ 
cipally  TE  modes  on  the  mode  sum  can  be  significant  at  frequencies  at  least  as 
low  as  20  kHz  for  the  case  of  westerly  propagation.  Azimuthal  anomalies 
included  drastic  polarization  changes  in  going  from  easterly  to  westerly 
paths.  In  the  case  of  transverse  propagation  at  the  magnetic  equator,  it  was 
shown  that  modes  which  are  pure  transverse  magnetic  (TM)  for  propagation  to 
the  east  may  be  pure  TE  for  propagation  to  the  west.  This  is  tantamount  to 
the  statement  that  modes  which  have  dominant  excitation  for  easterly  propaga¬ 
tion  may  have  vanishingly  small  excitation  for  westerly  propagation.  Azimu¬ 
thal  dependencies  are  characterized  by  rapid  variation  of  the  mode  constants 
in  the  neighborhood  of  north-south  or  south-north  propagation  both  for  equa¬ 
torial  and  midlatitudes.  These  variations  manifest  themselves  in  marked 
differences  in  mode  sums  for  azimuthal  changes  at  least  as  small  as  10°. 
Further,  it  was  shown  that  the  maximum  total  signal  attenuation  occurs  for 
north-south  propagation  rather  than  for  east-west  propagation. 

Although  Snyder  and  Pappert  (1969)  used  hypothetical  ionospheres,  their 
results  have  been  further  substantiated  by  Foley  et  al.  (1973)  using  "realis¬ 
tic  models  of  the  ionosphere."  Foley  et  al.  (1973)  essentially  duplicated  the 
analysis  of  Synder  and  Pappert  but  used  ionospheric  models  presented  in  the 
literature  by  Deeks  (1966)  and  Smith  et  al.  (1968).  Foley  et  al.  concluded 
that  mode  characteristics  are  dependent  on  the  ionosphere  model  assumed.  For 
example,  results  with  the  Smith  profile  were  similar  to  the  results  of  Snyder 
and  Pappert  while  under  certain  circumstances  results  with  the  Deeks  profile 
differed.  For  propagation  in  the  north-south  direction  with  the  Deeks  profile 
anomalies  were  found  in  the  behavior  of  modes  1  and  2.  It  would  appear  that 
between  16  and  20  kHz,  modes  1  and  2  interchange  their  identity.  Although 
unknown  to  Foley  et  al.  ,  this  identity  interchange  would  seem  to  be  just 
another  manifestation  of  the  modal  degeneration  reported  by  Snyder  (1968b). 
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2.  Observed  Propagation  Variations 


It  is  well  known  that  many  propagating  modes  may  be  supported  in  metallic 
waveguides  of  sufficiently  large  cross  section  compared  to  a  wavelength.  In 
such  a  guide  each  mode  has  a  different  phase  velocity,  producing  interference 
phenomena  as  the  various  modes  go  in  and  out  of  phase  at  varying  distances 
along  the  guide.  Furthermore,  changes  in  the  cross  section  of  the  guide  are 
known  to  produce  higher  order  modes  which  result  in  a  modified  modal  interfer¬ 
ence  pattern. 

The  earth-ionosphere  waveguide  is  four  to  six  wavelengths  in  height  at 
typical  VLF  wavelengths.  Higher  order  modes  can  be  excited  in  the  VLF  wave¬ 
guide  with  greater  amplitude  than  the  lowest  order  mode  even  thou<’’  the  lowest 
order  mode  suffers  less  attenuation.  Thus  modal  interference  as  a  function  of 
distance  from  a  VLF  transmitter  is  to  be  expected  just  as  for  a  metallic 
guide.  Such  is  indeed  the  case.  For  example,  VLF  field  amplitudes  were 
recorded  by  Rhoads  a  .d  Garner  (1967)  aboard  an  airplane  flying  between  Hawaii 
and  the  west  coast  of  the  United  States.  Comparison  of  their  daytime  measured 
amplitudes  with  results  computed  using  the  daytime  mode  constants  of  Wait  and 
Spies  (1964)  showed  almost  ideal  agreement.  Rhoads  and  Garner  (1967)  showed 
that  under  day  conditions  the  effect  of  the  earth’s  magnetic  field  on  the  VLF 
modal  constants  can  be  ignored,  at  least  for  distances  less  than  4000  km. 
They  also  confirmed  that  at  distances  up  to  at  least  3000  km  higher  order 
modes  must  be  considered.  The  good  agreement  obtained  for  daytime  data  could 
not  be  obtained  for  night  data. 

Early  investigations  of  VLF  propagation  to  great  distances  were  confined 
to  amplitude  measurements.  Because  of  the  distances  involved,  the  ground  wave 
could  not  be  used  as  a  reference,  and  oscillators  then  available  were  not 
sufficiently  stable  to  be  used  as  local  reference.  In  consequence,  it  was  not 
possible  to  study  the  phase  of  the  VLF  signals.  With  the  development  of 
atomic  frequency  standards  in  the  1960s,  investigation  of  the  phase  variabil¬ 
ity  of  VLF  signals  became  possible.  The  availability  of  such  stable  oscilla¬ 
tors,  combined  with  the  growing  interest  in  the  use  of  the  phase  of  VLF  sig¬ 
nals  for  navigation,  prompted  several  investigations  of  long  distance  propaga¬ 
tion  of  VLF  radiowaves. 
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Recall  that  one  of  the  characteristics  of  modal  excitation  and  attenua¬ 
tion  is  the  possible  change  in  modal  dominance  at  some  distance  from  a  trans¬ 
mitter.  At  the  point  where  mode  dominance  changes,  the  two  modes  are  of  equal 
amplitude.  If  at  this  point  the  two  modes  are  in  phase  opposition,  then  the 
resultant  mode  interference  will  result  in  a  deep  amplitude  minimum  accom¬ 
panied  by  local  instability  in  the  dependence  of  phase  on  distance.  For 
example,  a  minor  change  in  conditions,  such  as  a  small  change  in  height  of  the 
ionosphere,  results  in  a  stepwise  change  in  phase  by  as  much  as  2^.  Observa¬ 
tions  of  VLF  signal  phase  during  daytime  at  16  kHz  were  made  by  Burgess  (1967) 
aboard  an  inflight  aircraft.  A  significant  result  of  the  measurements  was  the 
detection  of  an  abrupt  (with  respect  to  distance)  change  in  phase  at  about 
3000  km  from  the  transmitter,  which  is  in  good  agreement  with  results  using 
the  modal  constants  from  Wait  and  Spies  (1964)  for  a  daytime  ionosphere. 

Another  manifestation  of  multimode  propagation  is  found  in  the  phenomena 
of  sunrise  and  sunset  fading.  This  fading  is  characterized  by  periodic  and 
repeatable  variations  in  amplitude  and  phase  as  the  dawn/dusk  terminator  moves 
along  a  VLF  propagation  path.  This  fading,  typically  most  pronounced  at  sun¬ 
rise  on  easterly  paths  and  at  higher  frequencies,  was  first  studied  by 
Yokoyama  and  Tanimura  (1933)  who  attempted,  without  success,  to  explain  their 
observations  in  terms  of  a  ray-optical  propagation  model.  Another  attempt, 
also  not  successful,  to  explain  the  fading  using  ray-optical  concepts  was 
presented  by  Rieker  (1963).  It  remained  for  Crombie  (1964)  to  explain  the 
phenomena  using  waveguide  mode  concepts.  Of  all  possible  combinations  of 
conditions  -  sunrise  and  sunset,  position  of  terminator  on  propagation  path, 
relative  positions  of  transmitter  and  receiver,  etc.  -  only  two  fundamental 
situations  occur  within  the  confines  of  Crombie 's  explanation. 

The  first  case  depicts  the  field  incident  on  the  terminator  from  the 
transmitter  as  determined  by  two  modes  and  the  field  at  the  receiver  deter¬ 
mined  by  only  one  irwde.  This  situation  would  be  characteristic  of  a  trans¬ 
mitter  on  the  night  side  of  the  terminator  and  the  receiver  on  the  day  side 
sufficiently  removed  to  avoid  multimode  propagation.  The  received  signal  then 
depends  on  the  amplitude  and  phase  relationships  of  the  modes  incident  on  the 
terminator  and  on  the  efficiency  with  which  these  modes  are  converted  by  the 
terminator  to  the  single  mode  which  reaches  the  receiver.  Because  the  modal 


relationships  are  determined  by  the  distance  of  the  terminator  from  the  trans¬ 
mitter,  as  the  terminator  moves  along  the  propagation  path  periodic  variations 
will  occur  in  the  received  signal.  The  periods  will  be  determined  by  the 
difference  in  the  phase  velocities  of  the  two  modes  entering  the  terminator, 
and  the  variations  must  occur  simultaneously  for  all  points  in  the  portion  of 
the  propagation  path  beyond  the  terminator  region.  The  fading  would  be  most 
intense  when  the  two  entering  modes  have  the  most  pronounced  interference, 
which  is  usually  when  the  terminator  is  near  the  transmitter. 

The  second  case  depicts  the  incident  field  as  a  single  mode.  The 
presence  of  the  terminator  leads  to  multimode  generation  in  the  region  past 
the  terminator  which  in  turn  leads  to  modal  interference.  In  this  case,  the 
interference  pattern  would  appear  to  be  attached  to  the  terminator  and  to  move 
synchronously  with  it.  Signal  minima  are  then  to  be  expected  at  fixed  dis¬ 
tances  from  the  terminator,  and  the  distance  between  fades  will  be  determined 
by  the  difference  between  the  phase  velocities  of  the  two  modes  leaving  the 
terminator  region.  The  deepest  signal  fades  would  be  expected  as  the 
terminator  approaches  the  receiver.  This  case  corresponds  to  propagation  from 
the  sunlit  side  of  the  terminator  into  the  night  side.  Note  that  in  both 
cases,  the  fade  spacing  would  be  determined  by  the  differences  between  the 
phase  velocities  of  two  nighttime  modes. 

A  thorough  experimental  examination  of  the  Crombie  model  was  conducted  by 
Walker  (1965)  who  recorded  18  kHz  transmissions  from  NBA,  Balboa,  Panama  while 
making  repeated  shipboard  crossings  of  the  Atlantic  in  the  equatorial  region. 
During  sunrise  (propagation  from  night  to  day)  signal  fades  were  observed 
simultaneously  at  all  points  in  the  day  portion  of  the  path  -  in  agreement 
with  the  first  case  discussed  above.  Walker  determined  the  change  in  termin¬ 
ator  position  between  successive  fades  and  from  this  calculated  the  required 
differences  in  the  phase  velocities  of  the  two  interfering  modes.  These 
results  were  in  good  agreement  with  theoretical  results  for  nighttime 
ionospheres  -  further  substantiating  the  Crombie  model.  Observations  were 
also  made  during  sunset  (propagation  from  day  to  night)  and  results  were  in 
excellent  agreement  with  the  Crombie  model. 
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3.  The  Equatorial  Anomal; 


The  observations  just  discussed  were,  for  the  most  part,  made  at  middle 
or  low  latitudes.  As  the  previous  discussion  indicates,  such  observations 
have  been  adequately  explained  using  quite  simple  VLF  waveguide  mode  theory. 
However,  some  observations  made  of  VLF  transmissions  over  long  transequatorial 
paths  have  apparently  indicated  an  anomalous  effect  associated  with  nighttime 
or  transitional  transequatorial  VLF  propagation. 

One  of  the  earliest  reports  on  transequatorial  VLF  propagation  to  show  an 
apparent  anomaly  was  by  Chilton  et  al.  (1964).  They  made  simultaneous  obser¬ 
vations  of  NBA  transmissions  (Canal  Zone,  Panama)  at  18  kHz  for  paths  to  both 
the  northern  and  southern  hemispheres.  Both  paths  were  of  similar  length  but 
the  southern  hemisphere  path  crossed  the  magnetic  equator.  Although  they  did 
not  look  at  transition  fading,  they  did  observe  an  anomalous  difference  in  the 
diurnal  changes  of  signal  amplitude  and  phase.  Chilton  et  al.,  suggested 
their  observations  resulted  from  a  difference  in  ionization  profile  due  to 
latitudinal  dependence  of  cosmic  rays.  More  recently,  Chilton  and  Crary 
(1971)  have  suggested  that  the  latitude-dependent  ionization  source  might 
arise  from  the  recently  discovered  x-ray  stars. 

Araki  et  al.  (1969)  recorded  dual  frequency  (15.5  and  22.3  kHz)  VLF  sig¬ 
nals  in  Japan  from  station  NWC  in  Australia.  They  reported  sunrise  fading  on 
22.3  kHz  but  not  on  15.5  kHz  and  further  noticed  that  during  fading  at  22.3 
kHz  the  first  phase  change  was  an  increase  followed  by  a  phase  decrease.  The 
nighttime  field  strength  was  somewhat  less  at  15.5  kHz  and  somewhat  higher  at 
22.3  kHz  than  during  the  day.  Araki  et  al.  (1969)  concluded  that  their 
results  suggest  the  existence  of  a  frequency-dependent  anomaly  of  the 
nighttime  field  relative  to  the  daytime  field  along  their  transequatorial 
propagation  path. 

Bickel  et  al.  (1970)  conducted  an  experiment  to  investigate  the  appar¬ 
ently  anomalous  results  of  Snyder  and  Pappert  (1969),  who  predicted  maximum 
signal  attenuation  in  a  southerly  rather  than  westerly  direction.  Measure¬ 
ments  of  23.4  kHz  signals  from  station  NPM  (Luaualei,  Hawaii)  were  made  aboard 
an  airplane,  flying  radials  along  the  propagation  paths  to  Seattle,  Ontario 
(California),  Samoa,  and  Wake  Island.  Bickel  et  al,  confirmed  the  Snyder  and 


1  3 


tfnnauaiNNi 


Pappert  (1969)  results  and  concluded  the  increased  southerly  attenuation  is  an 
effect  of  the  geomagnetic  field.  An  unexplained  effect,  however,  was 
observed.  For  the  radial  to  Samoa,  rapid  signal  amplitude  variations  with 
distance  were  observed  in  the  vicinity  of  the  geomagnetic  equator.  Although 
equipment  problems  may  be  suspected,  simultaneous  recordings  of  other  Hawaii- 
based  transmissions,  made  using  completely  separate  equipments,  displayed  sim¬ 
ilar  results  (Bickel,  private  communication,  1970). 

The  first  report  of  an  equatorial  anomaly  associated  with  transition 
fading  was  made  by  Lynn  (1967).  He  reported  on  VLF  transmissions  from  NLK 
(Jim  Creek,  Washington)  at  18.6  kHz  to  Smithfield,  South  Australia.  This 
propagation  path  is  approximately  13,500  km  long  and  the  receiver  to 
geomagnetic  equator  distance  along  the  path  is  approximately  5700  km.  The 
direction  of  propagation  is  essentially  southwesterly.  The  interference  dis¬ 
tance  reported  by  Lynn  for  sunrise  transition  fading  is  approximately  2000  km 
for  the  terminator  located  in  midlatitudes,  that  is,  in  excess  of  ±20°  from 
the  geomagnetic  equator.  Whenever  sunrise  transition  fading  was  observed 
while  the  terminator  was  within  ±20°  of  the  geomagnetic  equator,  the  interfer¬ 
ence  distance  increased  to  as  much  as  3700  km,  and  when  averaged  over  the 
anomaly  had  a  value  of  2900  km.  Lynn  (1967)  concluded  from  his  observations 
that  the  change  in  interference  distance  resulted  from  a  change  in  the  differ¬ 
ence  of  phase  velocity  for  two  modes  as  well  as  from  a  change  in  the  relative 
phases  of  the  appropriate  mode  conversion  coefficients.  The  question  of  a 
possible  cause  for  these  changes  was  left  unanswered  by  Lynn.  It  is  inter¬ 
esting  to  note  that  Lynn's  conclusions  are  dependent  on  his  selection  of  a  VLF 
propagation  model  as  well  as  his  selection  of  a  mode  conversion  model.  Lynn 
(1967)  assumed  an  isotropic  propagation  model  wherein  any  changes  in  propaga¬ 
tion  parameters  necessitates  changes  in  ionospheric  conditions. 

Kaiser  (1968)  has  examined  VLF  signals  from  NLK  at  18.6  and  24  kHz  and 
WWVL  at  20  kHz  received  at  Lower  Hutt,  New  Zealand.  These  propagation  paths 
are  similar  to  the  path  considered  by  Lynn  (1967)  and  Kaiser's  results  are 
much  the  same  as  Lynn's.  An  anomalously  large  interference  distance  (in 
excess  of  3000  km)  was  observed  for  sunrise  transition  fading  when  the  termin¬ 
ator  was  within  about  ±20°  of  the  magnetic  equator.  Kaiser  also  examined 
transition  fading  data  previously  published  elsewhere  and  concluded  that 


anomalously  large  values  of  the  sunrise  interference  distance  is  characteris¬ 
tic  of  generally  east  to  west  VLF  propagation  below  about  30°  magnetic  lati¬ 
tude.  A  further  conclusion  by  Kaiser  (1968)  is  that  any  anomaly  in  sunset 
transition  fading  for  such  paths  is  relatively  minor.  Kaiser  suggests  that  a 
possible  cause  of  the  anomaly  might  be  a  larger  and  sharper  day-night  wave¬ 
guide  transition  near  the  geomagnetic  equator  as  compared  to  higher  latitudes. 
He  does  not  however  suggest  a  possible  cause  for  the  latitudinal  dependence  of 
the  day-night  waveguide  transition. 

Additional  observations  of  this  transition  fading  anomaly  have  been 
reported  by  Lynn  (1969,  1970).  These  observations  further  substantiate  the 
general  characteristics  of  the  anomaly.  A  theoretical  interpretation  of  the 
transequatorial  anomaly  was  presented  by  Lynn  (1970).  His  interpretation  is 
based  only  on  assumed  latitudinal  variations  in  modal  phase  velocities  and 
relative  phases  of  the  mode  conversion  coefficients.  He  does  not  suggest  a 
possible  cause  for  the  variations. 

Meara  (1973)  has  extended  Lynn's  (1970)  analysis  to  include  determi nat ion 
of  the  changes  in  phase  velocities  within  the  equatorial  anomaly  region.  He 
concludes  that  the  phase  velocity  of  the  first  mode  is  essentially  unaltered 
when  propagating  through  the  equatorial  region,  whereas  that  of  the  second 
mode  is  reduced.  It  should  be  noted  that  Meara' s  results  are  entirely 
dependent  on  the  assumed  propagation  and  mode  conversion  models.  Once  again, 
no  explanation  of  a  possible  cause  for  the  variations  in  phase  velocities  is 
given. 

A  further  application  of  Lynn's  (1970)  analysis  to  explain  some  anomalous 
diurnal  changes  in  transequatorial  VLF  propagation  has  been  made  by  Arak  i 
(1973).  The  analysis  follows  essentially  that  previously  discussed.  Araki 
does  go  a  step  further  than  Lynn  (1970)  and  Meara  (1973)  however,  in  that  he 
claims  the  validity  of  the  isotropic  propagation  model  and  asserts  that  any 
change  in  propagation  parameters  results  from  a  modification  of  the  nighttime 
equatorial  ionosphere.  Lynn  (1978)  has  pointed  out  that  the  propagation  paths 
of  Chilton  et  al.  (1964)  and  Araki  (1973)  have  a  west-east  component  rather 
than  an  east-west  component,  as  do  the  other  paths  discussed  here.  He  further 
claims  that  only  a  minor  azimuthal  variation  in  relative  phase  velocities  of 
interfering  nighttime  modes  can  account  for  the  observations  of  Chilton  et  al. 


(1964)  and  Araki  (1973).  Lynn  (1978)  thus  concludes  that  the  observations  of 
Chilton  and  Crary  (1971)  cannot  be  seen  as  direct  evidence  for  the  control  of 
night  VLF  reflection  heights  by  stellar  x-ray  sources.  However,  Svennesson 
and  Wester lund  (1979),  using  waveguide  mode  theory  combined  with  nighttime 
profiles  given  by  ionospheric  theory,  claim  results  that  support  the  reports 
of  x-ray  star  effects  on  VLF  radio  signals. 

B.  Description  of  Research 

This  research  is  concerned  with  a  theoretical  and  observational  charac¬ 
terization  of  VLF  radiowave  propagation  in  equatorial  latitudes.  It  is  neces¬ 
sary  to  be  able  to  model  purely  propagational  effects  for  this  region  in  order 
to  interpret  observational  results  in  terms  of  changes  in  the  geophysical 
environment.  Examples  of  some  specific  questions  to  be  addressed  are: 

Can  the  anomalously  large  interference  distance  observed 
for  generally  east  to  west  transequatorial  paths  be  explained 
by  extending  conventional  VLF  propagation  models?  For  exam¬ 
ple,  can  the  anomaly  be  explained  simply  by  including  geomag¬ 
netic  field  effects  on  VLF  waveguide  mode  parameters  (as  sug¬ 
gested  by  Lynn)?  Can  the  equatorial  magnetic  field  per  se 
cause  enhanced  mode  conversion,  or  must  a  modified  terminator 
(as  suggested  by  Kaiser)  or  a  modified  nighttime  ionosphere 
(as  suggested  by  Araki)  be  employed  near  the  geomagnetic 
equator?  If  a  modified  equatorial  terminator  or  nighttime 
ionosphere  is  necessary,  what  causes  the  modification?  If 
geomagnetic  field  effects  on  the  modal  constants  is  the 
answer,  what  is  the  cause  for  the  effect? 

Trans-terminator  propagation  of  VLF  radiowaves  in  equatorial  latitudes 
for  generally  east  to  west  paths  is  characterized  by  the  anomaly  in  interfer¬ 
ence  distance.  In  order  to  address  the  question  of  a  cause  for  this  anomaly 
(e.g.,  geomagnetic  field  effects  on  waveguide  modal  parameters  vs.  a  modified 
equatorial  terminator),  data  are  needed  for  additional  propagation  directions. 
If  the  anomaly  is  caused  by  a  modified  equatorial  terminator,  it  seems  reason¬ 
able  to  expect  a  similar  effect  on  a  propagation  path  essentially  reciprocal 
to  the  paths  already  discussed.  If  it  is  caused  by  strictly  geomagnetic  field 
effects,  a  reciprocal  p^th  should  exhibit  dissimilar  results  or  perhaps  no 
anomaly  at  all. 


A  VLF  propagation  path  nearly  reciprocal  to  those  employed  by  Lynn  and 
Kaiser  is  available.  We  shall  analyze  VLF  signals  recorded  in  San  Diego  from 
station  NWC,  Australia.  This  provides  the  needed  information  for  a  reciprocal 

path. 

Relevant  theoretical  modeling  of  VLF  radiowave  propagation  will  be  per¬ 
formed.  This  will  include  consideration  of  inhomogeneities  both  in  the  ver¬ 
tical  direction  as  well  as  along  the  propagation  path.  Anisotropic  effects 
due  to  the  geomagnetic  field  will  be  considered. 
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II.  Theory 


A.  Horizontally  Homogeneous  Guides 

1 .  Coordinate  Transformation 

Formulation  of  waveguide  mode  theory  for  VLF  propagation  in  the  earth- 
ionosphere  waveguide,  homogeneous  in  and  transverse  to  the  direction  of 
propagation,  is  presented  in  this  section.  The  earth-ionosphere  guide  is 
modeled  as  cylindrical,  with  effects  of  curvature  transverse  to  the  direction 
of  propagation  accounted  for  through  use  of  an  appropriate  "spreading  term"  as 
described  in  Chapter  I.  Also  as  mentioned  in  Chapter  I,  the  formulation 
employed  in  this  study  is  based  on  the  pioneering  work  of  Budden  (1961, 
1962).  Because  Budden's  formulation  is  essentially  planar,  a  coordinate 
transformation  is  employed  to  convert  between  cylindrical  geometry  and  planar 
geometry.  The  transformation,  introduced  by  Richter  (1966)  with  respect  to 
earth  flattening  techniques,  and  used  by  Pappert  (1968)  and  Abbas  et  al. 
(1975)  in  VLF  propagation  studies  is  briefly  outlined. 

Application  of  the  conformal  coordinate  transformation  of  the  form  p  = 
aez^a;  4>  =  x/a  results  in  circles  (p  =  constant)  in  the  P ,  ■+>  plane  becoming 
straight  lines  (z  =  constant)  in  the  x,  z  plane.  If  a  represents  the  radius  of 
the  earth,  the  transformation  maps  the  surface  of  the  earth  (P  =  a)  to  the 
line  z  =  0.  The  transformation  further  maps  radii  ( ■}>  =  constant)  from  the  p ,  $ 
plane  to  lines  x  =  constant  in  the  x,z  plane.  A  full  circular  arc  (>  =  2t)  at 
P=0  is  mapped  so  that  the  length  of  x  is  the  circumference  of  the  earth. 

The  transformation  used  here  differs  from  that  employed  by  Richter 
(1966),  Pappert  (1968)  and  Abbas  et  al.  (1975)  in  that  they  maintained  some 
selected  altitude  (H)  as  invariant  in  both  coordinate  frames.  However,  the 
two  transformations  differ  only  in  quantities  of  second  order  in  smallness. 
As  the  invariant  altitude  (H)  approaches  zero,  the  two  transformations  become 
identical. 

Maxwell's  curl  equations  with  time  variation  of  the  form  elfl)t  for  a 
medium  characterized  by  a  dielectric  tensor  may  be  written  as 
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V  x  E  =  -ik  H 
~~  ~ 

V  x  H  =  ik  £  •  E 

~  ~  o*»  ~ 


where  kQ  is  the  wave  number  of  free  space  (kQ  =  u)/c)  and  £  is  the  dielectric 
tensor.  In  equations  (1),  the  variable  H  is  actually  the  magnetic  intensity 
multiplied  by  the  free  space  impedance.  Assuming  3/3y  =  0,  equation  (1)  in 
cylindrical  coordinates  may  be  written 
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Applying  the  coordinate  transformation  to  equation  (2)  yields 
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By  making  the  identifications 
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Ex  =  p/a  V  Ey  =  Ey;  Ez  =  p/a  EP; 
Hx  =  p/a  H<J>;  Hy  =  Hy'  Hz  =  p/a  V 


(4) 


equations  (3)  become  identical  to  what  would  be  obtained  in  a  rectangular 
coordinate  frame  for  a  medium  characterized  by  permeability  tensor  and 
dielectric  tensor  of  the  form 
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Some  additional  insight  is  obtained  by  considering  the  transformed  equa¬ 
tions  in  the  free  space  region  between  the  earth  and  ionosphere.  Making  the 
usual  substitutions  into  equations  (3c)  and  (3d)  from  equations  (3a)  and  (3b), 
along  with  the  free  space  assumption,  yields 
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Thus,  the  effect  of  the  transformation  is  to  replace  the  original  cylindrical 
free  space  region  with  a  planar  region  filled  with  a  medium  which  has  refrac¬ 
tive  index  equal  to  ez/a  for  TE  polarized  waves  or  which  has  permeability 
equal  to  e2Z//a  for  TM  polarized  waves. 
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2 .  The  Mode  Equation 


The  prototype  waveguide  model  considered  here  consists  of  two  plates  of 
infinite  extent,  parallel  to  the  x-y  plane  and  located  at  z  =  0  and  z  =  h. 
The  interior  of  the  guide  is  assumed  to  be  free  space.  This  prototype  wave¬ 
guide  corresponds  to  the  "fictitious  free  space  region"  to  be  later  discussed 
in  Section  4  where  it  will  be  used  to  determine  modal  excitation  factors  and 
polarization  coupling.  The  electromagnetic  field  is  assumed  to  be  independent 
of  the  y-coordinate  and  propagation  is  assumed  to  be  in  the  x-direction.  Any 
given  waveguide  mode  will  be  considered  as  composed  of  two  crossing  plane 
waves,  upgoing  and  downgoing,  with  wave  normals  at  the  angles  ±9  to  the  x- 
axis . 

The  nature  of  the  boundary  planes  will  be  specified  in  terms  of  reflec¬ 
tion  coefficients.  These  coefficients  are  defined  in  terms  of  the  ratios  of 
specific  components  of  the  upgoing  and  downgoing  waves  and  are  therefore 
dependent  on  the  polarization  of  the  waves.  A  convenient  way  of  describing 
the  polarization  is  in  terms  of  the  linear  components  of  the  complex  electric 
vector  which  are  perpendicular  and  parallel  to  the  plane  of  incidence  (the  x-z 
plane).  The  notation  11(1)  will  indicate  that  the  electric  vector  is  parallel 
(perpendicular)  to  the  plane  of  incidence. 

In  general,  the  upper  boundary  must  be  considered  anisotropic.  Thus, 
there  can  be  a  change  in  polarization  upon  reflection.  The  reflecting 
properties  of  the  upper  boundary  are  given  by  reflection  coefficients  having 
the  matrix  form 


R  = 
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where  the  first  subscript  describes  the  polarization  of  the  incident  (upgoing) 
wave  and  the  second  subscript  denotes  the  polarization  of  the  reflected  wave. 
The  properties  of  the  lower  boundary  are  assumed  isotropic  and  specified  by  a 
diagonal  reflection  matrix  of  the  form 
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By  considering  the  upgoing  wave  as  the  primary  or  incident  wave,  the 
downgoing  wave  can  be  considered  as  resulting  from  reflection  at  the  upper 
boundary.  This  wave  will  in  turn  be  reflected  at  the  lower  boundary  giving 
rise  to  another  upgoing  wave.  The  requirement  that  the  second  upgoing  wave  be 
the  same  as  the  original  yields  the  mode  condition.  This  process  is  shown  in 
the  diagram  below. 


The  primary  wave  is  represented  by  the  column  matrix  uq.  The  second  upgoing 
wave  ( )  is  therefore  given  by 

u  =  R  R  u  exp [-2ikhsin0] 

ss  1  S3  as  S3  Q  L  J 

where  the  exponential  term  accounts  for  phase  change  in  the  double  passage 
across  the  guide.  The  self-consistency  requirement  for  the  waveguide  mode  is 
thus 


det  [g  R  exp(-2ikhsin0)  -  1]  =  0  (10) 

where  J  is  the  unit  diagonal  matrix. 

3 .  Point  Sources  and  Line  Sources 

The  source  is  assumed  to  be  a  vertical  Hertzian  dipole.  Budden  (1961, 
1962)  has  shown  that  a  vertical  dipole  is  equivalent  to  an  infinite  set  of 
line  quadrupoles,  where  the  line  quadrupoles  are  located  parallel  to  the  x-y 
plane  and  inclined  at  (complex)  angles  to  the  y-axis.  This  is  easily  seer,  by 
considering  the  Hertz  vector  for  a  dipole  at  the  origin.  If  the  dipole  is 
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parallel  to  the  z-axis,  the  Hertz  vector  has  only  a  z-component  and  is  propor¬ 
tional  to  exp(-ikr)/r.  For  a  dipole  of  unit  strength,  the  Hertz  vector  is 
given  by  (Sudden,  1962) 
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where 

oi  -  a  ret  an  tz.\x  cos('j)  +  y  shi(j))  I, 

,  ,  ,  ,  a  z  1/2 

S  -  t x  cos(ui)  +  y  sin(j)  i  +  z 


Thus  s  is  the  distance  from  a  point  ;x,  y,  z)  to  a  line  in  the  x-y  plane 
passing  through  the  origin  at  an  angl-t  to  the  y-axis  and  a  is  the  angle 

between  the  x-y  plane  and  t.he  line  along  s.  is  a  Hankel  function  of  the 

second  kind  of  order  one.  The  term  cos  (  a)  H 1  ^2  *  ( ks )  would  result  from  a  line 
quadrupole  source  in  the  x-y  plane  at  an  angle  j  to  the  y-axis.  Consider  a 
line  quadrupole  source  as  composed  of  two  line  dipole  sources  at  x  -  V;>  $x  and 
x  =  -  V2  <Sx  in  the  x-y  plane.  The  line  dipoles  consist  of  closely  spaced 

dipoles  oriented  parallel  or  antiparallel  to  the  z-axis  and  distributed 
parallel  to  the  y-axis.  The  Hertz  vector  for  the  line  dipoles  is  proportional 
to  HQ^^(kp)  where  p  =  (  x2  +  z2 )  ^2  .  Letting  6x  tend  to  zero  and  letting  the 

strengths  (M,  -M)  of  the  line  dipoles  increase  in  such  a  way  that  the  product 

( 2  ) 

m6x  remains  finite  will  result  in  a  term  proportional  to  3h  (kp)/3x,  whicli 

o 

reduces  to  a  term  of  the  form  cos ( 0 ) H ^ 2 ) ( kp ) ,  where  cosO=x/p.  From  the  defi¬ 
nitions  of  s  and  a  in  (12),  it  is  clear  that  the  integral  in  (11)  results  from 
an  infinite  set  of  line  quadrupoles  in  the  x-y  plane,  each  at  a  different 

angle  w  to  the  y-axis.  Because  the  path  of  integration  in  ( 1 1 )  is  not  unique, 
there  is  an  infinite  number  of  ways  in  which  the  angles  u>  may  be  chosen. 


Consider  a  guide  bounded  by  perfectly  conducting  plates  with  a  line 

source  parallel  to  the  y-axis.  This  source  generates  waveguide  modes  and  each 

consists  of  crossing  plane  waves  with  normals  in  the  x-y  plane  at  angles  ±0m 

to  the  x-axis.  For  a  line  source  at  an  angle  w  to  the  y-axis,  the  component 

plane  wave  has  the  form  exp[-ik([x  cos(u>)  +  y  sin(w)]  C  +  zS  1  I  hr  r 

m  mJ  J  wnere  - 
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cos{0m)  and  Sm  -  sin(0m).  A  point  source  within  the  guide  is  new  represented 
by  an  infinite  set  of  line  sources  such  that  the  total  field  of  the  component 
wave  is  proportional  to 


Tt/2+i“> 

-n/2-i® 


exp[ -ikf  [x  cos(oj)  +  y  sin(w)]C  +  zS  )  1  dw. 

m  m  /  J 


(13) 


Making  the  substitutions  p 


<x< 


y2)  V2 


and  cos(u) 


x/p  in  (13)  gives 
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The  integration  in  (14)  yields  a  Hankel  function  of  the  second  kind  of  order 
zero  (Sommerfeld,  1949).  Thus,  the  total  field  of  the  component  wave  is 
proportional  to 


(2) 


H0  '(kpcm)  exp(-ikzSm). 


The  appropriate  Hankel  function  exchange  for  the  exponential  term  when  using  a 
line  source  is  therefore 


exp 


(-ikxc  )  + 
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4 .  Modal  Excitation  Factors  and  Polarization  Coupling 

The  excitation  of  waveguide  modes  is  conveniently  treated  in  terms  of 
line  sources.  The  technique  consists  of  selecting  one  of  the  possible  line 
sources  for  the  integrand  in  (11),  finding  the  mode  excitation  for  this  line 
source  and  then  effecting  a  change  from  an  exponential  term  to  the  appropriate 
Hankel  function.  Following  the  method  of  Budden  (1961,  1962),  the  source  is  a 
line  of  quadrupoles  at  height  d  (near  the  ground)  parallel  to  the  y-axis.  A 
fictitious  free  space  region  is  introduced  extending  from  z  =  1  to  z  =  u  such 


that  1  <  d  <  u  (the  line  source  is  assumed  within  the  free  space  region).  The 

region  above  z  =  u  is  characterized  by  the  reflection  coefficient  matrix  R 

»u 

while  the  region  below  z  =  1  is  characterized  by  a  reflection  coefficient 
matrix  g^.  Assuming  that  the  line  source  has  unit  strength,  the  y-component 
of  the  magnetic  field  in  infinite  free  space  can  be  written 

k  3 

O  t  TT  +  100  ? 

Hy  =  47TE-  J0-i»  exp[iko<xC+|z-d|S)]C  d0,  (16) 

1  o 

where  C  =  cos(9),  S  =  sin(O)  and  the  limits  of  integration  are  selected  to 
ensure  convergence  for  all  x  (Sommerfeld,  1949).  The  primary  field  integrand 
in  equation  (16)  represents  a  plane  wave  which  undergoes  reflections  at  z-u 
and  z=l,  and  these  reflected  waves  must  be  added  to  the  primary  field  to 
obtain  the  total  field.  Four  combinations  are  possible,  two  from  waves 
initially  upgoing  or  downgoing  at  the  source  and  two  from  waves  reaching  the 
reception  point  upgoing  or  downgoing. 

An  initially  upgoing  wave  is  converted  to  a  downgoing  wave  by  the 
reflection  process  denoted  by  the  operator  (g  g^)^g  and  to  an  upgoing  wave  by 
the  operator  (§j§u)  •  An  initially  downgoing  wave  is  converted  to  an  upgoing 
wave  by  the  operator  and  to  a  downgoing  by  the  operator  (g  g^)N. 

The  N  denotes  N  reflection  processes  of  the  the  type  described  by  the  bracket¬ 
ed  term.  Phase  terms  must  also  be  associated  with  each  reflection  process  to 
account  for  the  transversals  to-and-fro  between  z=l  and  z=u.  The  subscripts  u 
and  1  denote  reflection  coefficients  evaluated  at  z^u  and  z=1,  respectively. 

The  total  field  is  obtained  by  combining  the  primary  and  all  reflected 
fields  and  letting  the  dimensions  of  the  free  space  region  tend  to  zero.  This 
requires  taking  the  limits  1+u  and  z+d  in  such  a  way  that  1  <  d  <  u  is  main¬ 
tained.  The  total  field  may  be  written  as 
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o  N=  1 


The  reflection  coefficients  are  evaluated  at  z=d  and  the  subscripts  l,u  are  no 
longer  needed.  In  equation  (17),  the  first  summation  starts  with  N=0  in  order 
to  include  the  primary,  nonreflected  wave,  and  P  is  a  column  matrix  whose 
second  element  is  zero  corresponding  to  a  vertical  electric  dipole  source. 


The  order  of  summation  and  integration  is  now  reversed  so  that  each 
integrand  will  contain  a  geometric  series.  It  can  be  shown  (Budden,  1961, 
1962)  that  the  series  converges  for  certain  simplified  cases,  and  convergence 
is  assumed  here.  Summation  of  the  geometrical  series  gives 
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Consequently  the  result  of  summing  the  series  in  equation  (17)  is 
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If  the  sign  of  0  is  reversed,  the  role  of  incident  and  reflected  wave  is 
exchanged  and  each  reflection  element  is  replaced  by  its  inverse.  The 
integrand  simply  changes  sign  upon  reversal  of  the  sign  of  9.  The  integral 
may  be  evaluated  using  the  theory  of  residues  by  deforming  the  contour  to  run 
from  n/2  +  i°°  to  -ir/2  -  i°°  symmetrically  through  the  origin.  Because  the 

integrand  is  antisymmetric  about  9=0,  the  integral  along  the  new  contour  is 
zero.  In  changing  the  contour,  however,  singularities  are  crossed,  the  most 
important  of  which  are  poles  of  the  factor  | 1  -  |  g  )  .  The  residue  contri¬ 

bution  of  a  simple  pole  at  9  is 
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where 


F  =  det  [£  -  §  g] 

X  =  lim  {F( 1  -  I  §]_1 
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and  the  subscript,  n,  signifies  that  the  quantities  are  evaluated  at  9  -  Qn< 
The  pole  condition  F  =  0  is  identical  to  the  waveguide  mode  equation  (10)  in 
the  limit  h  *  0 .  Thus,  the  solution  of  the  pole  condition  yields  the  wave¬ 
guide  modes  and  the  residues  at  the  poles  give  the  excitation  factors. 
Performing  the  matrix  multiplications  in  equation  (22),  using  the  Hankel 
function  substitution  (equation  (15))  and  employing  the  geometrical  spreading 
factor  yields 
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To  be  compatible  with  the  mode  sum  of  Chapter  I  (equation  (1))  it  is  necessary 
to  change  the  reference  direction  for  0  from  the  x-axis  to  the  z-axis.  Doing 
so  simply  requires  the  interchange  of  sinQ  and  cosG. 

Equation  (24)  requires  the  source  and  reception  point  to  be  at  the  same 
height,  which  is  also  the  height  at  which  the  reflection  elements  are  to  be 
determined.  If  the  dependence  of  Hyn  on  height  z  is  given  by  the  function 
Gn(z),  then  an  obvious  extension  of  equation  (24)  to  include  arbitrary  source 
and  reception  height  involves  multiplication  by  the  function 


G  (z  ) 
n  T 


G  (d) 


G  (d)  ' 


where  zT  is  the  transmitter  height,  zR  is  the  receiver  height  and  d  is  the 
evaluation  height  for  the  reflection  elements,  assumed  to  be  below  the 
ionosp.ere. 

Whenever  the  ionosphere  is  assumed  to  be  isotropic  (no  geomagnetic  field) 
the  field  components  can  be  divided  into  two  sets.  These  are  the  well  known 
TE  ( Hy  =  0)  and  TM  (Ey  =  0)  modes.  When  anisotropy  is  included,  however,  a 
given  mode  in  general  involves  both  Hy  and  Ey.  To  determine  the  coupling 

between  the  TM  and  TE  components  of  the  eigenfunctions  again  we  introduce  a 

thin,  fictitious  region  of  free  space  within  the  transformed  planar  waveguide. 
The  region  is  assumed  to  be  centered  at  some  altitude  z=d  which  is  below  the 
ionosphere  so  that  the  region  z<d  can  be  considered  isotropic.  within  this 
region,  the  fields  of  the  TM  and  TE  components  can  be  considered  in  terms  of 

upgoing  and  downgoing  plane  waves.  Furthermore,  the  eigenfunctions  satisfy 

the  condition  that  a  particular  component,  say  the  upgoing,  must  be  returned 
to  its  original  value  after  being  reflected  from  the  two  boundaries  of  the 
thin  free  space  region.  Thus,  for  the  upgoing  component  e^  (the  component 
perpendicular  to  the  plane  of  propagation),  the  requirement  is 
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In  the  free  space  region 
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where  h^u  is  the  upgoing  component  of  Hy  and  e^11  is  the  upgoing  component  of 
Ey.  The  total  Hy  and  Ey  fields  within  the  free  space  gap  at  z-d  are  therefore 
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and  their  ratio  is 
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In  equation  (26)  the  reflection  coefficients  are  evaluated  at  z=d. 

5 .  Determination  of  Reflection  Coefficients 

Two  basic  methods  have  been  employed  for  the  numerical  calculation  of 
ionospheric  reflection  coefficients  for  horizontally  stratified  anisotropic 
media.  One  involves  the  numerical  integration  of  a  set  of  differential  equa¬ 
tions  describing  the  EM  fields  directly,  such  as  the  equations  of  Clemmow  and 
Heading  (1954).  Use  of  the  equations  of  clemmow  and  Heading  implies  consider- 
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ation  of  the  ionosphere  as  a  cold  magnetoplasma.  Instability  problems  are 
encountered  if  the  field  equations  are  directly  integrated.  This  may  be 
avoided  by  introducing  a  new  variable  for  integration,  which  does  not  have 
such  instabilities. 

Thus,  Budden  (1955)  considered  integration  of  the  linear  Clemmow  and 
Heading  equations,  pointing  out  the  numerical  instabilities,  and  then  intro¬ 
ducing  a  non-linear  differential  equation  for  the  reflection  coefficients.  It 
can  be  easily  shown  that  the  differential  equation  for  the  reflection  coeffi¬ 
cient  is  of  the  Riccati  type  for  many  simple  ionospheric  models.  As  pointed 
out  by  Ralston  and  Wilf  (1967),  the  Riccati  equation  exhibits  excellent  sta¬ 
bility  for  numerical  integration.  Budden  (1955)  pointed  out  that  the  new 
integration  not  only  did  not  have  numerical  instabilities  but  also  required 
less  computer  time  and  core  than  integration  of  the  wavefields  directly. 
Barron  and  Budden  (1959)  introduced  a  different  variable,  the  admittance 
matrix,  which  they  claimed  required  less  computer  time  and  core  than  the 
reflection  matrix.  However,  use  of  the  admittance  matrix  sometimes  leads  to 
numerical  difficulties.  A  disadvantage  of  using  either  the  reflection  or 
admittance  matrix  is  that  it  is  difficult  to  obtain  the  wavefields  directly 
from  the  functions  employed  in  the  integration. 

Pitteway  (1955)  introduced  a  numerical  way  to  overcome  the  instability. 
Pitteway's  techniaue  was  to  integrate  the  Clemmow-Heading  field  equations 
directly  and  to  use  Gram-Schmidt  or thogona 1 izat ion  at  various  steps  in  the 
integration  process.  An  advantage  of  this  scheme  is  that  wavefields  within 
the  ionosphere  are  made  available.  However,  this  is  achieved  at  the  price  of 
large  computer  core  requirements.  Alternative  schemes  using  Pitteway's  basic 
idea  have  been  developed.  Walsh  (1967)  and  Seliga  (1966),  for  example,  inte¬ 
grate  a  set  of  equations  for  what  they  term  "vacuum  modes"  (Poeverlein,  1967), 
applying  Gram-Schmidt  orthogona lization  to  avoid  numerical  instability. 

The  other  basic  method  involves  dividing  the  inhomogeneous  ionosphere 
into  thin  homogeneous  layers.  For  example,  Johler  and  Harper  (1962)  devised  a 
matrix  of  large  dimension  to  describe  a  multislab  ionosphere,  and  Price  (1964) 
simplified  the  matrix  to  a  series  of  4x4  matrices.  The  problem  of  numerical 
instability  was  not  addressed.  Wieder  (1968)  attempted  to  solve  the  numerical 
xnstability  problem  with  the  Price  4x4  matrix  by  performing  the  iterative 
multiplication  from  the  bottom  of  the  ionosphere  upward,  stopping  where 
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reflection  coefficients  seemed  to  stabilize.  Although  usually  effective,  this 
scheme  is  also  subject  to  numerical  problems.  Nagano  et  al.  (1975)  solved  the 
numerical  swamping  pri  ,lem  using  Gram-Scbmidt  orthogonalization  with  the 
iterative  matrix  multiplication  applied  to  a  layered  ionosphere. 

The  original  reflection  coefficient  matrix  integration  (Budden,  1955), 
using  standard  fourth  order  Runge-Kutta  techniques  ( Ralston  and  Wilf,  1967), 
has  been  adopted  for  this  study.  The  integration  is  started  at  some  heiubt 
above  which  negligible  reflection  occurs.  Error  control  and  adjustment  of  the 
integration  step  size  is  determined  by  comparing  second-order  Runge-Xutta 
results  with  fourth-order  results. 

Analytic  expressions  are  obtained  for  the  ground  reflection  coefficients. 
Consider  looking  downward  at  soma  height  (d)  below  the  ionosptiere  in  the 
transformed  planar  guide  of  Section  1.  in  order  to  be  compatible  with  the 
assumptions  of  Sections  2  and  4,  the  region  above  d  i.s  assumed  to  be  free 
space,  whereas  the  region  below  d  is  filled  with  the  "equivalent"  permeability 
or  dielectric  medium  as  described  in  Section  1.  The  fields  in  the  region  be¬ 
tween  z-d  and  z=0  satisfy  equations  (7a, b)  while  the  fields  above  d  and  within 
the  earth  have  simple  exponential  forms. 

Recognizing  that  all  altitudes  of  interest  are  small  compared  to  the 

radius  of  the  earth  (a),  the  exponential  term  in  equation  (7)  is  expanded  to 
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first  order  in  z/a.  A  variation  in  the  x-direction  of  the  form  e  °  is 
assumed.  Note  that  later  the  variable  S  will  be  identified  as  sin(6)  where 
the  angle  0  is  referenced  to  the  vertical  rather  than  the  horizontal  as  in 
previous  sections.  With  the  above  assumptions,  equation  (7b)  for  TM  waves 
becomes 


H"  +  k2  (C2  +  ctz)  H  =  0,  UD 
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where  C  =  1-S  ,  a  =  2/a  and  the  double  prime  denotes  second  differentiation 

with  respect  to  z. 

With  the  change  of  variable 
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the  equation  becomes 
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where  now  the  double  prime  denotes  differentiation  with  respect  to  t.  Equa¬ 
tion  (29)  is  Stokes'  equation.  Thus,  H  is  linearly  expressible  in  terms  of 
solutions  to  Stokes'  equation.  In  this  study,  we  shall  use  the  modified 
Hankel  functions  of  order  one  third  (h^,h?),  as  defined  by  the  Computation 
Laboratory  at  Cambridge,  Massachusetts  (1945).  Using  a  similar  procedure  for 
TE  polarization,  we  shall  arrive  at  similar  results  except  that  the  trans- 
formed  planar  guide  will  contain  a  medium  with  refractive  index  given  by  e  '  . 

The  waves  in  the  three  regions  are  given  by 
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Region  II: 
0  <  z  <  d 


Region  III: 
z  <  0 
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where  ^ R jj  is  the  TM  reflection  coefficient  being  sought,  and  are  modi¬ 
fied  Hankel  functions  of  order  one  third,  N  is  the  refractive  index  of  the 
ground  and  T  =  ( N2  -  S2)1^2,  taken  with  negative  imaginary  part.  Continuity 
requirements  between  the  regions  leads  to  the  results 
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Following  a  similar  procedure,  results  for  the  TE  polarization  are  given  by 
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B.  Horizontally  Inhomogeneous  Guides 

If  significant  variations  occur  along  the  direction  of  propagation,  which 
is  often  the  case  in  the  earth- ionosphere  waveguide,  then  homogeneous  guide 
formulations  require  modification.  For  such  cases,  there  are  two  techniques 
for  mode  summation  that  are  commonly  used  at  VLF.  One  employs  a  WKB  approach 
and  the  other  uses  mode  conversion. 

In  the  WKB  technique,  an  eigenmode  is  assumed  to  be  uniquely  and  inde¬ 
pendently  identifiable  anywhere  along  a  propagation  path.  Furthermore,  each 
mode  is  assumed  to  depend  only  on  the  local  characteristics  of  the  waveguide 
and  to  propagate  independently  of  the  existence  of  any  other  modes.  The  WKB 
formulation  generalizes  the  expression  for  ttie  mode  summation  in  equation  (1), 
Chapter  I,  to  the  form 
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Note  that  the  only  differences  arise  from  the  use  of  the  geometric  mean  of  the 

T  R 

excitation  factors  at  transmitter  (A  )  and  from  receiver  (A  )  locations  and 

n  n 

from  the  integration  of  the  complex  propagation  factor  along  the  path.  The 
WKB  formulation  was  apparently  introduced  simply  as  an  "obvious  extension"  of 
the  homogeneous  mode  sum  by  Wait  (1964),  who  further  noted  that  the  form  is 
similar  to  what  would  be  expected  if  a  WKB  formulation  were  applied,  although 
no  actual  derivation  was  presented.  The  use  of  the  name  WKB  mode  sum  has 
persisted. 

In  the  mode  conversion  method,  eigenmodes  are  assumed  to  be  uniquely  and 
independently  identified  only  in  a  local  sense.  Even  though  the  eigen 
characteristics  of  a  mode  depends  only  on  local  characteristics  of  the 
waveguide,  the  propagation  of  each  mode  depends  not  only  on  the  existence  of 
other  modes  but  also  on  the  past  and  future  propagation  history  of  itself  and 
the  other  modes. 

The  approach  using  mode  conversion  does  not  lend  itself  conveniently  to 
consideration  of  transition  fading  on  very  long  VLF  propagation  paths.  The 
WKB  method,  although  easily  applied,  is  unable  to  handle  the  rapid  change  in 
guide  characteristics  that  occur  at  the  terminator  (Pappert  and  Snyder, 
1972).  A  difficulty  encountered  when  using  mode  conversion  arises  from  the 
necessity  of  knowing  the  future  propagation  history  of  each  mode.  Making 
appropriate  assumptions  to  relieve  this  requirement  can  lead  to  some  simplifi¬ 
cation,  but  the  resulting  formulation  is  still  not  convenient  for  very  long 
paths.  A  further  difficulty  is  that  conventional  mode  conversion  techniques 
require  a  knowledge  of  all  the  vector  components  of  the  electromagnetic  field 
throughout  all  space. 

Numerous  methods  have  been  presented  for  circumventing  this  latter 
problem.  Wait  ( 1968a, b)  for  example,  formulated  mode  conversion  using  isotro¬ 
pic  impedance  boundary  conditions  with  the  result  that  fields  were  required 
only  for  the  region  between  the  upper  (ionospheric)  impedance  boundary  and  the 
lower  (ground)  impedance  boundary.  Such  boundary  conditions  are  certainly  not 
appropriate  for  the  nighttime  ionosphere  and  thus  not  for  consideration  of 
transitions  from  night  to  day.  Smith  (1974)  extended  Wait's  (1968a, b)  results 
to  include  boundary  anisotropy  by  using  impedance  matrix  boundary  conditions. 
His  analysis  however,  was  for  a  planar,  empty  waveguide  and  was  directed 
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towards  considering  ground  effects.  For  the  ionosphere  Smith  (1974)  assumed 
the  quasi-longitudinal  approximation  (given  by  Budden,  1961)  for  a  homogeneous 
cold  plasma,  and  this  is  not  suitable  for  the  study  at  hand. 


Galejs  (1971)  presented  a  modified  mode  conversion  formulation  for  an 
abrupt  change  between  daytime  and  nighttime  anisotropic  ionospheres.  Conven¬ 
tional  mode  conversion  theory  usually  represents  the  EM  field  as  an  expansion 
in  terms  of  a  suitability  chosen  complete  set  of  orthogonal  functions.  The 
Galejs  formulation  employs  an  incomplete  set  of  non-or thogona 1  functions.  The 
formulation  now  to  be  employed  is  similar  to  that  of  Galejs  (1971). 


1 .  The  Operator  Equations 

The  operator  equations  presented  in  this  section  are  based  on  Maxwell's 
equations  in  two  (rectangular)  dimensions  and  are  patterned  after  the  formula¬ 
tions  of  Pappert  and  Smith  (1972).  Furthermore,  material  media  are  assumed  to 
be  characterized  by  constitutive  parameters  (permittivity  and  permeability) 
that  can  be  inhomogeneous  in  one  dimension.  This  is  quite  compatible  with  the 
usual  procedure  for  handling  the  earth-ionosphere  waveguide  by  describing  the 
electromagnetic  properties  of  the  ionosphere  with  the  aid  of  a  vertically 
inhomogeneous  dielectric  tensor. 


We  start  with  equation  (3),  use  the  substitutions  in  equation  (4),  and 
eliminate  the  Ex  and  Hx  components.  This  leads  to  the  differential  operator 
equation 


fe  e  = 


ik  ax' 
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(33) 


where  e  is  a  four  element  column  vector  given  by 


(34) 


and  ^  is  a  4x4  matrix  differential  operator.  Equation  (33)  is  cast  as  an 

eigenvalue  equation  by  using  the  functional  form  exp( - i kQsnx ) ,  so 

that  L  e  =  S  e  .  These  e  form  the  set  of  expansion  functions  used  in 
a  — n  n— n  ~n 

conventional  mode  conversion  theory. 
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To  find  an  operator  adjoint  to  L  a  procedure  discussed  by  Pappert  and 
Smith  (1972)  is  used.  First,  let  the  inner  product  be  defined  by 

<a,b>  =  /  a*  •  b  dz  (35) 

where  the  tilde  over  a  denotes  transpose  and  the  asterisk  denotes  conju¬ 
gation.  Next,  the  operator  adjoint  to  L  is  defined  by  the  relationship 

+ 

<L  w,e>  2  <w,Le>  (36) 

which  is  just  a  prescription  for  integration  by  parts  to  find  the  form 
of  L  as  well  as  the  necessary  boundary  conditions  on  the  adjoint  functions 
(w)  in  order  that  the  integrated  contributions  vanish.  Pappert  and  Smith 
(1972)  showed  that  the  adjoint  functions  can  be  obtained  from  direct  functions 
by  considering  a  second  waveguide  which  is  related  to  the  original  one  in  a 
simple  way.  In  the  original  guide,  the  z-axis  is  taken  as  vertical,  and 
propagation  is  assumed  in  the  xz-plane.  The  x-direct ion-cosine  of  the 
geomagnetic  field  is  1  in  the  original  guide.  The  adjoint  guide  is  obtained 
by  replacing  1  with  -1.  By  examining  symmetry  properties  of  the  elements  of 
the  ionospheric  permittivity  tensor,  Pappert  and  Smith  (1972)  showed  that  the 
adjoint  eigenfunctions  are  given  in  terms  of  eigenfunctions  of  the  adjoint 
waveguide  according  to  the  relationship 
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Further,  from  the  definition  of  adjoint  given  by  equation  (36)  the  usual 
biorthogonality  relations  hold.  Thus 


<L  w  ,e  > 
=>  ~n  ~m 


=  X*  <w  ,e  > 

n  ''n  ~m 

-  <w  ,L  e  > 
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which  shows  that  <w  ,e  >  =  0  as  well  as  X  -  S* ,  for  nondeqenerate  eigenfunc- 
~n  ~m  n  n 

tions . 

2 •  Mode  Conversion 

Fundamental  to  the  mode  conversion  formulation  employed  here  is  the 
notion  that  any  variation  in  the  VLF  propagation  conditions  along  a  path  can 
be  modeled  as  a  series  of  discrete  changes  between  homogeneous  sections  along 
the  path.  With  this  in  mind,  it  is  seen  that  the  fundamental  analytical 
requirement  is  the  determination  of  mode  conversion  coefficients  due  to  a 
single  abrupt  change  in  waveguide  characteristics. 

Following  classical  procedures,  consider  a  mode  of  unit  amplitude  inci¬ 
dent  on  such  a  discontinuity  in  a  waveguide.  This  mode  is  described  in  terms 

>  j 

of  its  height-gain  function  e^.  (z),  where  the  superscript  I  denotes  the  inci¬ 
dent  region,  the  subscript  j  denotes  the  jth  mode,  the  overarrow  indicates  a 
forward  propagating  mode,  and  the  undertilde  indicates  a  vector  quantity. 
Note  that  the  terms  height-gain  function  and  eigenfunction  are  used  here 
interchangeably  even  though  more  traditional  usage  limits  the  term  height-gain 
to  individual  components  of  the  vector  eigenf unction.  The  discontinuity  in 
the  guide  generates  both  reflected,  or  backward  propagating,  modes  within  the 
incident  region  (I)  and  transmitted,  or  forward  propagating,  modes  within  the 
transmission  region  (II).  Backward  propagating  modes  might  also  exist  in 
region  II  due  to  another  discontinuity;  these  will  therefore  be  included.  The 
continuity  requirements  lead  to  the  equation 
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(38) 


which  is  independent  of  the  position  of  the  discontinuity  between  regions  I 

and  II.  In  equation  (38),  the  coefficient  rmj  is  a  coefficient  describing  the 

conversion  of  the  incident  mode  (j)  into  a  backward  traveling  mode  (m)  in 

region  I  (eI(z)),  and  T_  ^  is  a  coefficient  describing  conversion  of  the 

n  j 

incident  mode  (j)  into  the  forward  traveling  mode  (n)  of  region  II  (e^3(z)). 

~n 

The  R^j  are  the  coefficients  appropriate  to  possible  backward  traveling  modes 
in  region  II  (e**(z))> 


37 
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Recall  that  a  major  problem  with  implementation  of  mode  conversion 
concepts  arises  from-  the  need  to  know  the  future  propagation  history.  In 
terms  of  equation  (38),  this  concept  is  represented  by  the  coefficients  R^ . 
By  making  the  assumption  that  appreciable  reflection  does  not  occur,  so  that 
there  can  be  no  backward  traveling  modes,  the  R^  may  be  set  equal  to  zero. 
The  rmj  are  then  also  all  zero  and  equation  (38)  reduces  to 


+  I  „  +TT 

e  .  (z)  =  I  T  .  eAi(z)  .  (39) 
~j  n  nj  ~n 

The  assumption  of  negligible  reflection  is  not  unique  to  this  study.  It  is 
commonly  made  and  has  been  demonstrated  to  be  adequate  in  numerous  cases 
(Wait,  1968a, b;  Galejs,  1971;  Pappert  and  Snyder,  1972;  Smith,  1974).  In 
equation  (39),  there  are  assumed  to  be  N  modes  in  region  II  so  that  there  are 
N  conversion  coefficients. 

Conventional  mode  conversion  techniques  requirf  multiplication  of 
equation  (39)  by  the  transpose  conjugate  of  the  adjoint  function  for  the  mth 
forward  traveling  mode  in  region  II,  followed  by  integration  over  all  space. 
Doing  so,  and  recalling  the  biorthogonality  relationships,  we  obtain  the  well- 
known  result 
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more  conveniently  written  as 
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Solving  for  the  conversion  coefficient  yields  che  apparently  simple  result 
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(41  ) 


The  simplicity  of  equation  (40)  depends  on  knowing  the  height-gain  func¬ 
tions  throughout  all  space.  Such  a  requirement  is  the  main  limitation  to  the 
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application  of  mode  conversion  techniques  to  very  long  VLF  paths.  Relaxing 
the  integration  limits  to  some  finite  value  prevents  use  of  the  biorthogona¬ 
lity  relationships  and  yields,  instead  of  equation  (40),  the  result 
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where  there 

are  N  equations 

and 
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In  equation  (42 
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where  the  overarrow  notation  is  dropped  because  only  forward  propagating  modes 
are  considered. 

For  application  to  the  VLF  earth- ionosphere  problem,  the  above  procedure 
may  be  justififed  as  follows.  First,  the  various  eigenvalues  and  eigenfunc¬ 
tions  are  determined  using  the  complete  geometry,  (-“  <  z  <  +»).  Second,  the 
eigenfunctions  can  be  expected  to  decay  rapidly  within  the  earth  and  within 
the  ionosphere.  Below  the  ionosphere  but  within  the  waveguide,  the  vector 
eigenfunctions  degenerate  into  the  isotropic  TM  and  TE  components  with  a 
prescribed  coupling  between  the  two  for  a  given  mode.  By  taking  the  upper 
limit  •  1  the  integral  to  be  within  the  guide  and  below  the  ionosphere  it  is 
possible  to  write  the  various  In  ^  in  rms  of  known  functions;  namely  the 
modififed  Hankel  functions  discusse  1  in  Section  A. 5.  The  procedure  used  here 
is  actually  expansion  in  terms  of  "non-orthogona 1 ”  functions  and  is  certainly 
better  than  assuming  orthogonality  and  "throwing  away"  the  contributions  to 
the  integrals  from  the  ionosphere. 

The  jth  vector  eigenfunction  for  either  region  is  defined  to  be  (from 
equation  ( 34  )  ) 
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Because  the  upper  limit  on  the  integrals  in  equation  (43)  is  to  be  taken  below 
the  ionosphere,  equation  (44)  can  be  written  as 
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The  adjoint  vector  eigenfunction  is  given  by  equation  (37)  as 


h*  (-1)  V 

/  S*  1 

zm  1 

’  m 

-h*  (-1) 

-h* 

ym 

ym 

-e*  (-1) 

S*  1 

zm 

m 

e*  (-1)  / 

\  e* 

ym 

ym 

where  the  notation  indicates  that  elements  of  the  adjoint  eigenfunction  are 
obtained  by  reversing  the  x-axis  direction-cosine  of  the  earth's  biasing  mag¬ 
netic  field.  The  first  integral  in  equation  (43)  thus  becomes 
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From  Section  A. 4  and  A. 5  it  is  seen  that  h  (-1)  =  h  but  that  e  (-1)  * 
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eym'  However,  from  equation  (36)  of  Section  A. 4 
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Because  -1 )  =  ^  and  (R^(-1)  =  (Budden,  1955,  1961b),  equation  (46) 

becomes 
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Denoting  this  ratio  as  A^  equation  (45)  becomes 
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The  integrals  in  equation  (48)  are  now  easily  evaluated  by  using  the  con¬ 
ventional  "multiply  and  subtract"  method.  Using  the  method  with  equation  (27) 
for  inodes  of  index  m  and  n,  which  may  be  in  the  same  region  or  differing 
regions,  leads  to  the  result 
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The  second  integral  in  equation  (48)  thus  becomes 
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while  the  first  integral  in  equation  (48)  is 
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In  equations  (49)  and  (50),  it  is  assumed  that  Cm11  *  for  all  m  and  j. 

Evaluation  of  the  second  integral  in  equation  (43)  proceeds  in  the  same 
manner.  The  only  difference  is  the  fact  that  all  field  components  in  this 
case  pertain  to  region  II  of  the  waveguide.  In  this  case  however,  the 
situation  m=j  results  in  division  by  zero  in  equations  (49)  and  (50)  and 
special  treatment  is  required.  For  the  case  m=j  in  region  II,  the  use  of  L' 
Hospital's  rule  yields 
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where  the  superscript  "II"  is  omitted. 

The  individual  height  gains  can  be  written  from  equations  (30),  (31)  and 

(26)  as 
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Once  the  conversion  coefficients  are  obtained,  it  is  necessary  to  sum 
over  all  the  incident  n»des.  This  requires  evaluation  of  the  amplitudes  of 
the  incident  modes.  This  can  be  done  by  starting  at  the  transmitter,  where 
the  modal  excitation  factors  give  the  mode  amplitudes,  and  using  the  assump¬ 
tion  that  the  regions  between  the  discontinuities  are  homogeneous. 
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Notationally  this  can  be  written  as 


A1  = 


I  I 
P  a 


where 


represents  the  amplitudes  of  the  incident  modes  (A.)  at  the  end  of  the  region 
I  and 


represents  the  amplitudes  of  the  modes  at  the  start  of  region  I.  The 
factor  g1  is  a  diagonal  matrix  with  the  nth  diagonal  element  given  by 
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Note  that  Pn  represents  the  propagation  factor  appropriate  to  a  homogeneous 
section  cf  length  (Ap)1,  assumed  to  be  the  length  of  region  I.  Note  that  the 
amplitudes  of  the  modes  in  region  IX  just  after  conversion  can  be  conveniently 
written  as 
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where  J  is  a  rectangular  matrix  of  the  conversion  coefficients  from  region  I 
into  region  II. 

A  mode  conversion  formulation  has  l>een  presented  by  Pappert  and  Shockey 
(1974),  but  important  differences  exist  between  their  treatment  and  that 
employed  here.  This  occurs  not  only  in  the  numerical  procedures  employed  but 
also  in  the  definition  of  the  conversion  coefficients.  In  the  present  formu¬ 
lation  the  conversion  coefficient  is  computed  at  each  slab  interface  indepen¬ 
dently  of  any  previous  conditions.  In  the  Pappert  formulation,  the  conversion 
coefficients  are  defined  and  computed  in  such  a  way  that  the  effect  of  the 
propagation  factor  (equation  (55))  is  included.  Both  definitions  of  conver¬ 
sion  coefficients  are  reasonable,  although  their  numerical  values  can  be  quite 
different.  For  example,  first  consider  the  case  of  an  isolated  terminator 


specified  as  an  abrupt  transition.  The  two  formulations  will  give  the  same 
results  for  the  conversion  coefficients.  Next  suppose  that  the  terminator  has 
a  substantial  extent.  The  conversion  coefficients  then  differ  for  each  slab 
interface  after  the  first  because  the  Pappert  formulation  includes  the  accum¬ 
ulated  effects  of  propagation  through  the  transition  region,  whereas  the 
current  formulation  does  not. 


The  Pappert  definition  of  conversion  coefficients  is  useful  for  studying 
a  horizontal  inhomogeneity  existing  over  a  small  distance,  such  as  the  dawn- 
dusk  terminator.  A  similar  conversion  coefficient  can  be  derived  from  the 
mode  constants  and  the  mode  conversion  coefficients  presented  in  this  study. 
What  is  required  is  the  complex  amplitude  of  a  mode  at  the  end  of  such  a 
confined  inhomogeneity  resulting  from  a  single  mode  of  unit  amplitude  at  the 
start  of  the  inhomogeneity. 


Denote  by  region  1  the  region  prior  to  the  start  of  the  inhomogeneity. 

The  complex  amplitudes  of  the  modes  in  the  first  region  of  inhomogeneity 

(called  region  2)  due  to  a  mode  of  unit  amplitude  and  number  k  in  region  1  are 

2  1 

given  by  the  conversion  coefficients  T  (see  equation  (42)).  The  complex 

m,  k  £ 

amplitudes  of  the  modes  at  the  end  of  region  2  are  A  where 
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and  p^  is  the  propagation  factor  for  mode  m  in  region  2  (see  equation  (55)). 

The  complex  amplitude  of  mode  n  at  the  start  of  the  next  region  (region  3)  due 

3  2 

to  a  mode  of  unit  amplitude  in  region  2  is  T  '  ,  so  that  the  actual  complex 

n ,  m 

3  2 

amplitude  of  the  mode  n  in  region  3  is  given  by  the  product  of  T  ' ^ 
and  A2  Summation  over  all  the  modes  of  region  2  gives  the  total  complex 

amplitude  of  mode  n.  Thus 
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(57) 


This  result  can  be  generalized  to  any  inhomogeneity  along  the  path.  We  obtain 


ap  =  £  pP-1  P-' 

m,  k  n  m,  n  n  n,  k 
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The  formulations  presented  in  this  Chapter  have  been  implemented  in  a 
computer  program  written  in  FORTRAN.  Before  discussing  some  computed  results, 
however,  we  shall  discuss  observations  of  signals  received  in  San  Diego  from 
station  NWC  in  Australia. 


III.  EXPERIMENTAL  RESULTS 


The  majority  of  the  experimental  results  relating  to  the  transequator ia 1 
anomaly  discussed  previously  in  Chapter  I  pertain  to  propagation  paths  which 
are  essentially  southwesterly.  Further,  some  theoretical  results  presented  in 
Chapter  I  suggested  a  marked  difference  between  the  transequatorial  propaga¬ 
tion  to  the  southwest  as  compared  to  the  northeast.  Transequatorial,  trans¬ 
terminator  propagation  in  a  northeasterly  direction  is  thus  expected  to  be 
different  from  that  for  a  southwesterly  direction.  In  order  to  determine  the 
validity  of  this  expectation,  experimental  results  are  needed. 

A  seemingly  ideal  situation  exists  for  obtaining  such  results.  A  VLF 
communications  station  is  located  at  the  North  West  Cape,  Australia.  This 
station  is  currently  broadcasting  a  phase  stable  signal  at  22.3  kHz,  a  fre¬ 
quency  sufficiently  high  to  insure  the  occurrence  of  mode  conversion  at  the 
ionospheric  terminator.  A  receiving  site  located  at  San  Diego,  California 
further  presents  an  ideal  situation  in  that  the  propagation  path  from  North 
West  Cape  to  San  Diego  is  nearly  reciprocal  to  the  paths  employed  by  Lynn 
(1969).  NWC  transmissions  recorded  in  San  Diego  during  1975  are  discussed  in 
this  chapter. 

A.  Equipment  Used 

The  NWC  transmissions  were  received  at  San  Diego  using  essentially  "off 
the  shelf"  items,  commercially  available.  In  particular,  one  receiver  was  an 
EECO  Model  880  VLF  Receiver  and  two  were  TRACOR  Model  599  VLF  Tracking 
Receivers,  factory  modified  to  provide  amplitude.  Integration  time-constants 
for  all  receivers  were  from  one  to  three  minutes,  with  three  minutes  most 
jften  used.  All  receivers  obtained  phase  reference  from  a  cesium  clock  stan¬ 
dard. 

one  of  the  receivers,  the  EECO,  was  connected  most  of  the  time  to  a 
commercially  available  loop  antenna  and  part  of  the  time  to  a  whip  antenna. 
One  of  the  TRACORs  (no.  2)  was  always  connected  to  a  whip  antenna  and  the 
other  TRACOR  (no.  1)  was  connected  to  a  Textran  Model  611  Cardioid  Unit.  The 
cardioid  unit  was  used  to  "null  out"  any  possible  long  path  signal  from  the 
NWC  transmitter. 


Signal  levels,  in  dB  above  an  arbitrary  reference,  were  recorded  on 
Easterline  Angus  Recorders  (Model  S6015)  running  at  .75  inches  per  hour, 
resulting  in  a  time  resolution  of  about  two  minutes.  The  signal  phase  from 
the  EECO  was  recorded  using  a  Rustrak  Model  88  Recorder  whereas  Easterline 
Angus  Recorders  were  used  with  the  Tracors.  Two  receivers  were  always  tuned 
to  NWC,  TRACOR  no.  1  and  either  the  EECO  or  the  other  TRACOR.  The  redundancy 
of  receivers  tuned  to  NWC  was  to  provide  a  "back-up"  in  case  of  failure  in  one 
of  the  receivers,  a  not  infrequent  occurrence. 

The  receivers  were  calibrated  by  feeding  a  rf  signal  through  a  Hewlett 
Packard  Model  355D  Attenuator  into  the  receiver  antenna  input  from  an  EECO 
Model  883  Coherent  Frequency  Synthesizer.  The  chart  unit  deflections  for  a 
given  input  voltage  were  recorded  and  converted  to  an  equivalent  dB  level. 
All  signal  levels  are  thus  known  on  a  relative  but  not  an  absolute  scale. 

An  additional  calibration  procedure  was  followed  occasionally.  All 
receivers  were  tuned  to  the  same  frequency,  usually  22.3  kHz  for  NWC,  for  at 
least  one  full  diurnal  period.  Any  differences  were  noted  in  order  to  ascer¬ 
tain  if  the  various  antenna  systems  used  were  influenced  differently  by 
changing  local  factors.  Such  differences  were  never  found  to  be  significant 
except  during  rainy  days  when  the  whip  antenna  apparently  became  flooded. 

Conversion  of  data  from  analog  to  digital  form  was  accomplished  through 
use  of  a  Talos  graphics  tablet.  This  tablet  is  an  electromagnetic  sensing 
probe  and  plate  with  approximately  a  2  ft  by  2  ft  surface.  Changes  in  analog 
data  to  within  ±0.005  inches  can  be  digitized  with  the  graphics  tablet.  A 
"reading"  line  was  inscribed  on  the  analog  tracings  of  the  EA  recorders. 
Because  the  more  than  60  dB  dynamic  variation  in  received  signal  amplitude 
exceeded  the  full  scale  deflection  capabilities  of  the  recorders,  the 
anplitudes  of  deep  signal  fades  were  often  lost.  The  time  of  occurrence  of 
the  minima,  however,  could  be  accurately  determined. 

The  digitized  data  were  processed  by  computer  to  determine  times  of 
signal  minima.  Signal  amplitudes  were  also  processed,  but  they  have  not  been 
used  in  this  study  due  to  the  frequent  loss  of  minima  data.  For  preliminary 
processing  a  data  point  was  considered  to  correspond  to  a  signal  minimum  if  it 
was  more  than  0.25  dB  below  both  the  three  previous  and  three  following  data 
points.  Further  examination  of  the  preliminary  data  was  then  made  in  order  to 
eliminate  any  erroneous  data.  Data  obtained  during  periods  when  only  one 
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receiver  was  operating  properly  were  examined  with  particular  care.  Before 
discussing  the  final  processed  data  some  properties  of  the  propagation  paths 
and  the  ionospheric  sunrise/sunset  terminator  are  discussed. 


B.  The  Terminator  and  the  Propagation  Path 

The  ionospheric  sunrise/sunset  model  used  in  this  study  is  the  one  com¬ 
monly  employed  in  VLF  propagation  analysis.  Shown  in  figure  la*  is  a  pictor¬ 
ial  representation  of  the  terminator  model.  The  model  depicts  an  ozone  layer 

of  thickness  hn  which  is  assumed  to  be  a  perfect  absorber  of  relevant  solar 
u3 

radiation.  The  terminator  model  shows  that  ionospheric  sunrise  and  sunset  are 
given  in  terms  of  solar  zenith  angles  according  to  the  expression 

a  +  hQ 

3 

sinv  ,  -  - - -  i  =  day  or  night  , 

1  a  +  h  . 


where  a  is  the  radius  of  the  earth  and  h^  is  usually  assumed  to  be  70  km  for 
day  and  90  km  for  night.  Assuming  a  nominal  50  km  ozone  layer  leads  to  the 
results  xD  *  94°,  xN  “  96°  •  Typically,  only  the  value  xN  is  used. 

Another  important  property  of  the  transition  region  model  is  shown  in 
figure  1b.  In  that  figure  is  a  plan  view  of  the  transition  region,  with  an 
intersecting  propagation  path.  Note  that  the  length  of  the  propagation  path 
within  the  transition  is  given  essentially  as 


I. 


Lmin 


/cosB 


where  planar  geometry  is  assumed  to  be  adequate  for  the  distances  involved. 
If  the  above  values  of  xD  and  x^  are  used,  a  value  of  about  200  km  is  obtained 
for  Lm^n.  Even  though  the  value  of  Lm£n  depends  on  a  choice  of  solar  zenith 
angles  for  the  day  and  night  sides  of  the  transition  region,  note  that  the 
ratio  L/Lm^n  depends  only  on  &•  A  pictorial  representation  of  the  intersec¬ 
tion  of  the  propagation  path  with  the  transition  region  is  shown  in  figure  2. 


•Figures  appear  at  the  end  of  the  document. 
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The  period  during  which  a  particular  propagation  path  is  in  total  day, 
total  night  or  transition  depends  on  a  number  of  factors,  including  path 
orientation  and  location,  time  of  year,  etc.  Periods  of  night,  day  and 
transition  for  the  NWC  to  San  Diego  path  are  shown  in  figure  3  using  an 
ionospheric  terminator  corresponding  to  a  solar  zenith  angle  of  96°  .  Note 
that  the  propagation  path  is  predominantly  transitional  over  any  24  hour 
period  throughout  the  entire  year.  Maximum  sunrise  transition  period  is 
during  the  northern  hemisphere  summer  while  for  sunset  it  is  during  the 
winter. 

C.  Transition  Fading  Data 

1 .  Time  of  Signal  Minima 

The  times  of  minimum  signal  strength  are  shown  in  figure  4  for  NWC 
transmissions  at  22.3  kHz  as  received  at  San  Diego,  California  during  night  to 
day  transition  (sunrise).  The  relative  terminator  movement  between  transmit¬ 
ter  and  receiver  is  depicted  in  the  upper  right  insert.  In  this  case  the 
terminator  moves  from  receiver  location  to  transmitter  location  so  that  the 
first  fade  (in  time)  occurs  with  maximum  night  length  to  the  path.  The  last 
fade  occurs  with  minimum  night  length  to  the  path.  Seven  signal  minima  are 
observed.  The  data  of  figure  4  indicate  the  variability  of  the  times  of  sig¬ 
nal  minima  throughout  the  year.  Maximum  length  of  sunrise  fading  occurs  near 
midyear,  consistent  with  the  data  of  figure  3.  Daily  variability  is  also 
shown  by  figure  4.  The  times  of  signal  minima  fluctuate  by  as  much  as  5  to  10 
min  from  day  to  day.  The  effect  of  averaging  the  times  of  the  minima  over  a 
15  day  period  are  shown  in  figure  5. 

Times  for  sunset  signal  minima  are  shown  in  figure  6.  Also  shown  is  the 
relative  terminator  movement  between  transmitter  and  receiver.  In  this  case, 
the  first  fade  occurs  with  maximum  day  length  to  the  path.  Seven  signal  mini¬ 
ma  are  indicated  for  the  sunset  transition  as  for  the  sunrise.  The  variabil¬ 
ity  of  fading  times  throughout  the  year  again  are  consistent  with  the  data  of 
figure  3.  The  day  to  day  variability  of  the  sunset  data  is  noticeably  larger 
than  for  the  sunrise  data.  Times  of  sunset  siqnal  strength  minima  averaged 
over  15  days  are  shown  in  figure  7. 


2.  Terminator  Position  at  Maximum  Signal  Fade 


Associated  with  the  times  of  occurrence  of  minima  is  the  location  of  the 
terminator  at  those  times.  The  location  of  the  terminator  is  here  defined  to 
be  the  point  along  the  propagation  path  where  the  solar  zenith  angle  is  96°. 
A  more  precise  definition  would  be  the  point  on  the  propagation  path  where  the 
angle  from  the  sub-solar  point  is  96°.  The  transmitter/terminator  distance  at 
times  of  signal  minima  for  night-to-day  transition  (sunrise)  were  computed 
from  the  data  of  figure  4  and  are  shown  in  figure  8.  The  fade  numbering 
scheme,  based  on  time,  is  shown  at  the  right.  Also  shown  in  the  figure  are 
the  relative  positions  of  the  geomagnetic  equator  and  20°  magnetic  latitude 
(north  and  south).  Signal  minima  4  and  5  are  well  within  the  equatorial  region 
and  minima  3  and  6  are  marginally  so.  Thus,  if  the  equatorial  spacing  anomaly 
in  transition  fading  is  to  occur  for  this  path,  the  data  should  clearly  so 
indicate. 

The  distances  moved  by  the  sunrise  terminator  between  signal  minima, 
shown  in  figure  9,  are  nearly  constant  throughout  the  year,  with  the  possible 
exception  of  the  minima  1  and  2.  The  changes  in  these  minima  could  possibly 
be  due  to  the  proximity  of  the  transition  region  to  the  receiver  location, 
resulting  in  interference  effects  from  additional  modes.  Even  without  any 
smoothing  or  averaging,  the  data  contain  no  anomalies  due  to  the  presence  of 
the  terminator  in  the  equatorial  region. 

Transmitter/terminator  distances  for  times  of  signal  minima  for  day-to¬ 
night  transition  (sunset)  are  shown  in  figure  10.  There  are  significant 
differences  from  the  sunrise  case.  There  is,  first  of  all,  an  obvious 
seasonal  dependence.  This  dependence  does  seem  to  be  common  to  all  the 
fades.  A  perusal  of  the  original  strip  charts  showed  no  indications  of 
equipment  failure.  It  seems,  then,  that  there  is  a  significant  seasonal 
difference  between  the  sunrise  and  sunset  fading. 

There  is,  however,  an  even  more  significant  difference  between  the 
sunrise  and  sunset  fades.  This  is  indicated  in  fades  numbered  4,  5,  6  and 
7.  During  early  midyear,  there  seems  to  be  a  change  of  as  much  as  1500  km  in 
the  location  of  the  terminator  for  the  fades  4,  5,  h  and  7.  Note  that  these 
fades  are  essentially  within  the  equator i a 1  region,  thus  suggesting  the  possi¬ 
bility  of  an  equatorial  anomaly  for  the  sunset  transition  case. 


So 


An  explanation  of  the  apparent  anomaly  may  be  found  by  considering  the 
relative  length  of  the  transition  region  throughout  the  year  as  shown  in 
figure  11.  The  ratio  L/Lm^n  is  computed  as  1/cos3,  where  3  is  the  angle 
depicted  in  figure  1b.  Note  from  the  lower  figure  (for  sunset)  that  the 
effective  length  of  the  transition  region  increases  rapidly,  reaching  a 
maximum  near  midyear.  This  increase  in  the  effective  length  of  the  transition 
region  has  apparently  reduced  the  mode  conversion  effects  so  greatly  that  the 
"fades"  are  little  more  than  ripples  in  the  amplitude  data.  This  is  seen  in 
figure  12,  which  shows  the  amplitude  during  sunset  for  June  1,  2,  and  3,  1975. 
In  figure  12,  at  best,  four  transition  fades  can  be  identified. 

The  data  presented  here  have  demonstrated  no  geomagnetic  field  anomaly 
associated  with  the  northeasterly  propagation  direction.  Any  apparent 
anomalies  can  be  explained  in  terms  of  purely  geometrical  effects.  There 
remains  yet  the  unexplained  differences  in  seasonal  effects.  It  would  seem 
reasonable  to  search  for  an  explanation  of  the  seasonal  effects  in  terms  of 
the  seasonal  behavior  of  the  ionosphere,  a  subject  not  considered  in  th ' s 
study . 


IV.  NUMERICAL  MODELING  RESULTS 


Empirical  observations  discussed  in  Chapter  I  indicate  an  apparent 
anomaly  in  the  nighttime  VLF  propagation  across  the  geomagnetic  equator  for 
southwesterly  paths.  Limited  theoretical  results,  also  discussed  in  Chapter 

I,  suggested  such  anomalies  might  not  occur  for  northeasterly  propagation 
across  the  geomagnetic  equator.  VLF  propagation  data  taken  during  this  study 
confirmed  this  latter  point.  Using  the  theoretical  model  liscussed  in  Chapter 

II,  numerical  modeling  results  have  been  obtained  for  both  northeasterly  and 
southwesterly  transequatorial  VLF  propagation.  These  results  are  discussed  in 
this  chapter. 

A.  Review  of  VLF  Propagation  Environment 

A  crucial  input  to  the  theoretical  determination  of  VLF  propagation  is 
the  environmental  description.  Long  distance  propagation  of  VLF  radio  waves 
is  considered  herein  by  mode  theory  in  which  the  earth  and  the  lower  part  of 
the  ionosphere  are  considered  to  be  the  boundaries  of  a  lossy  waveguide.  For 
this  reason,  a  review  of  the  VLF  propagation  medium  requires  a  consideration 
of  the  electrical  properties  of  the  earth  and  the  lower  ionosphere. 

The  earth  is  usually  considered  to  be  a  homogeneous  isotropic  conductor. 
For  most  of  the  surface,  the  effective  ground  conductivity  ranges  from  about  a 
miilimho/m  for  typical  soils  to  about  5  mho/m  for  sea  water.  For  large  areas 
covered  with  perma-frost  or  ice,  such  as  the  Greenland  ice  cap,  the  conductiv¬ 
ity  can  be  as  low  as  10“'’  mho/m.  For  conductivity  this  low,  the  skin  depth 
can  be  so  large  that  consideration  must  be  given  to  the  lower  substrata 
(Westerlund,  1974). 

The  detailed  geometrical  characteristics  of  the  surface  are  not  very 
important  for  VLF  waves  because  the  principal  electric  vector  is  largely 
perpendicular  to  the  surface.  However,  if  the  electric  vector  has  an  appreci¬ 
able  horizontal  component,  these  characteristics  car.  be  quite  important 
(Galejs,  1972).  yuite  the  opposite,  however,  is  true  concerning  the  conduc¬ 
tivity  (Westerlund,  1974). 

The  upper  boundary  of  the  VLF  waveguide  is  considerably  more  complicated 
than  the  lower.  The  reflection  of  VLF  waves  from  the  upper  boundary  (the  D- 
region  and  lower  E-region)  is  generally  accepted  to  lie  from  about  70  km  for 
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daytime  and  about  90  km  for  nighttime.  The  reflecting  medium  is  a  slightly 
ionized  plasma,  rendered  anisotropic  by  the  geomagnetic  field.  Because  of  the 
variable  solar  control,  the  variations  in  electro-chemical  processes  and  the 
variations  in  thermal  and  mechanical  influences,  the  medium  is  further  ren¬ 
dered  highly  inhomogeneous. 

The  reflection  process  is  quite  a  complicated  process,  of  which  all 
aspects  have  still  not  been  fully  explained.  Both  the  density  gradients  and 
collision  frequencies  are  known  to  play  a  major  role  (Budden,  1961).  Reflec¬ 
tion  at  VLF  does  not  occur  at  the  altitude  where  the  operating  frequency  is 
suitably  related  to  the  plasma  frequency,  as  in  HF  ionospheric  propagation, 
but  is  controlled  to  a  significant  extent  by  the  ratio  of  the  plasma  to  the 
collision  frequency.  It  can  be  said,  however,  that  the  upper  boundary  of  the 
earth-ionosphere  waveguide  is  controlled  by  the  ionization  at  heights  below 
about  100  km. 

It  would  be  nice  to  be  able  to  report  that  the  D-  and  lower  E-regions  are 
well  understood.  This  is  not,  however,  the  case.  It  has  been  asserted,  that 
many  of  the  most  interesting  unsolved  problems  of  the  ionosphere  are  related 
to  the  D-region  (Bowhill,  1975).  In  a  recent  review  of  the  D-region,  Thomas 
(1974)  pointed  out,  that  although  there  have  been  advances  in  theoretical 
models  of  the  D-region  in  recent  years,  there  has  not  been  a  significant 
improvement  in  our  understanding  of  the  aeronomical  processes  operating. 

Ionization  below  100  km  is  produced  by  electromagnetic  and  corpuscular 
radiation  which  bombards  the  earth  from  a  variety  of  sources.  The  principal 
ionization  sources  for  the  quiet  D-region  are  considered  to  be  galactic  cosmic 
rays,  solar  x-rays  and  solar  H  Lya,  both  direct  and  scattered.  An  additional 
source  only  recently  recognized  and  not  yet  universally  accepted  involves 
precipitation  of  energetic  electrons. 

The  photoionization  of  nitric  oxide  by  solar  H  Lya  constitutes  the  major 
source  of  ionization  during  the  daytime  between  approximately  65  and  90  km 
altitude.  Soft  x-rays  (31-100  A)  dominate  above  90  km  and  galactic  cosmic 
rays  are  most  important  below  65  km.  At  nighttime,  ionization  is  usually 
assumed  to  be  produced  primarily  by  H  Lya  in  the  nightglow.  A  drizzle  of 
energetic  particles  is  continually  raining  down  onto  the  midlatitude  atmos¬ 
phere.  According  to  Potemra  and  Zmuda  (1970),  the?  electrons  are  the  domi¬ 
nant  nighttime  ionization  source  above  about  80  km.  Additional  sources  of 


nighttime  ionization  have  been  proposed,  including  cosmic  x-rays  and  solar  UV 
ionization  of  metastable  oxygen.  These  have,  however,  been  discounted  or 
demonstrated  to  be  of  only  secondary  importance. 


The  immediate  effect  of  these  ionization  sources  is  the  production  of 
free  electrons  and  positive  ions  of  atomic  or  molecular  size.  The  electrons 
can  combine  with  neutral  particles,  forming  negative  ions.  Positive  ions  and 
electrons,  as  well  as  positive  and  negative  ions,  can  recombine  forming  neu¬ 
trals.  Balance  equations  describing  these  processes  can  be  formulated  and 
solved.  However,  the  coefficients  involved  in  such  balance  equations  are 
usually  quite  uncertain.  Additional  chemical  and  physical  processes,  such  as 
the  presence  of  ions  of  larger  size  than  single  molecules,  further  complicate 
the  balance,  resulting  in  the  frequent  use  of  "nominal  working  models"  of  the 
ionosphere.  Such  models  range  from  a  homogeneous  conductor  sharply  bounded  at 
an  assumed  reflection  height,  to  an  "electron  only"  (positive  ions  assumed  to 
be  infinitely  massive)  exponentially  varying  ionosphere,  to  a  "lumped  param¬ 
eter"  model  involving  electrons  together  with  positive  and  negative  ions  of  a 
single  mass. 

Ionospheric  profiles  of  the  lumped  parameter  type  were  recently  presented 
by  Rawer  et  al.  (1978)  for  electrons.  These  profiles,  designated  by  the  name 
IRI-78,  are  obtained  from  computer  models  that  attempted  to  represent  averages 
of  measured  ionospheric  parameters.  A  set  of  exponentially  varying  electron 
density  profiles  were  recently  recommended  for  VLF  propagation  studies  by 
Mot  citt  (1977).  These  profiles  were  obtained  by  fitting  multifrequency  long 
path  VLF  propagation  data  to  theoretically  computed  data  on  a  tr ial-and-error 
basis.  An  extensive  comparison  of  the  two  sets  of  profiles  was  made  in  order 
to  ascertain  which  to  employ  in  this  study.  Typical  comparative  results  are 
displayed  in  figures  13  and  14. 


In  figure  13,  a  mode  sum  computed  using  the  IRI-78  nighttime  profile  and 
magnetic  conditions  corresponding  to  propagation  from  station  NPM  in  Hawaii  to 
Ontario,  California,  is  presented  along  with  VLF  amplitude  data  recorded 


aboard  an  airplane  flying  on  such  a  path.  The  computed  mode  sum  differs 
significantly  from  the  data.  The  corresponding  "best  fit"  exponential  elec¬ 
tron  density  profile  is  given  by  Morf itt  (1977)  as  one  witti  reference  height 
{ h ’  )  of  85.5  km  and  gradient  (6)  of  0.5  km-1.  The  B,h'  notation  is  adopted 

2/v  =  2.5  x  105  exp( 3( z-h' ) ) , 

from  Wait  and  Spies  (1964)  who  assumed  that  Up 


where  z  is  the  height  above  ground  in  kilometers,  the  angular  plasma 

frequency,  and  v  is  the  electron  collision  frequency,  given  by  v  =  1.816  x 
1011  exp  [-.15z]  .  Thus,  a  change  in  3  or  h'  denotes  a  change  in  electron 
density  only.  In  figure  14,  similar  results  for  a  daytime  case  are  presented. 
The  mode  sum  computed  using  the  IRI-78  is  in  better  agreement  in  this  case 
than  for  the  nighttime  case  but  is  nonetheless  still  not  in  very  good  agree¬ 
ment.  The  "best  fit"  exponential  profiles  in  the  daytime  case  are  3  =  0.3,  h’ 
=  75  for  3  February.  Based  on  numerous  comparisons  such  as  these,  exponential 
electron  density  profiles  are  employed  in  this  study.  The  nighttime  profile 
selected  has  3  =  0 .5  km-1  and  h*  =  86  km  while  the  daytime  has  3  =  0 .3  km-1 
and  h'  =  74  km. 


B.  Sample  Calculation  Comparisons 

Height  gains  computed  using  the  program  based  on  the  theory  of  Chapter  II 
have  been  compared  with  those  used  in  the  work  of  Pappert  and  Snyder  (1972). 
Three  electron  density  profiles  are  considered:  3  =  0.5  km-1,  h  =  86  km  and 

70  km  and  3  =  0.3  km-^,  h  =  70  km.  The  geomagnetic  field  is  taken  to  be 
horizontal,  corresponding  to  the  magnetic  equator,  and  the  propagation 
direction  is  taken  to  be  directly  west  to  east  so  that  there  is  no  coupling 
between  the  TM  and  TE  eigenfunctions.  Only  the  TM  eigenfunctions  are 
considered  here.  The  earth  is  assumed  to  be  sea  water  and  a  frequency  of  18.6 
kHz  is  used.  These  conditions  correspond  to  those  used  by  Pappert  and  Snyder 
(1972) . 

Shown  in  figure  15  are  the  magnitudes  of  three  eigenfunctions  for  the  3  = 
0.5  km  h  =  70  km  profile.  The  right  hand  part  of  the  figure  shows  the  h^ 
components  whereas  the  left  hand  part  shows  the  ex  component.  The  .lashed 
curves  for  altitudes  greater  than  about  65  km  are  derived  from  Pappert  and 
Snyder  (1972)  and  demonstrate  the  decay  of  the  eigenfunctions  within  the 
ionosphere.  Below  65  km  there  is  complete  agreement  between  the  results 
obtained  here  and  those  of  Pappert  and  Snyder,  thereby  verifying  the  results 
obtained  in  this  study. 

The  profile  here  is  relatively  abrupt  in  that  the  region  of  the  ionos¬ 
phere  from  about  65  to  70  km  is  the  most  important.  The  region  below  65  km  is 
essentially  free  space  whereas  above  70  km  the  waves  are  largely  evanescent 
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and  decay  by  about  two  powers  of  ten  from  their  nominal  value  within  the 
guiding  region. 

Results  for  a  less  abrupt  profile  (3  =  0.3  tan-1,  h  =  70  tan)  are  shown  in 
figure  16.  From  the  results  of  Pappert  and  Snyder  concerning  phase  velocity 
(their  figures  3  and  4),  it  follows  that  the  effective  vertical  width  of  the 
waveguide  for  this  profile  is  about  the  same  as  for  the  3  =  0 .5  km-1  profile. 
However,  as  can  be  seen  from  figure  16,  the  guide  is  essentially  free  space 
only  below  about  55  or  60  km. 

Mode  conversion  coefficients  have  been  computed  using  the  formulation 
discussed  in  Chapter  II.  For  this  purpose,  region  1  was  assumed  to  have  the  3 
=  0.5  km1;  h  =  86  tan  profile  and  region  2  was  assumed  to  have  the  3  =  0.3 
km  h  =  70  km  profile.  This  situation  corresponds  to  one  reported  by 

Pappert  and  Snyder.  The  coefficients  for  converting  incident  modes  1  and  2 
into  transmitted  mode  1  are  given  in  table  1  using  upper  limits  of  integration 
from  55  to  70  km.  The  other  conversion  coefficients  were  insignificant, 
consistent  with  the  results  of  Pappert  and  Snyder.  Best  agreement  is  obtained 
for  an  upper  limit  of  integration  equal  to  60  km.  The  results  computed  here 
agree  with  the  more  exact  results  of  Pappert  and  Snyder  to  within  about  30%  to 
40%.  The  electromagnetic  fields  computed  in  region  2  by  the  methods  presented 
here  differ  f rom  those  of  Pappert  and  Snyder  by  only  a  few  dB. 

Table  1.  Magnitude  of  the  Mode  Conversion  Coefficients 
Conversion  from  Incident  M  to  Transmitted  N 

Mode  Number  Integration  Limit 


N 

M 

70 

65 

60 

55 

Exact  Results 

1 

1 

2.75 

2.02 

1 .94 

2.28 

1  .47 

1 

2 

1.10 

0  .76 

0  .72 

0  .89 

0  .53 

C.  Modeling  Results 

The  numerical  modeling  results  will  be  presented  in  two  parts.  The  fiist 
part  addresses  the  northeasterly  t ransequator ia 1  path.  The  propagation  path 
from  station  NWC  (22.3  kHz),  Australia  to  San  Diego,  California  is  used  for 
this  case.  The  second  part  addresses  the  southwesterly  path  from  station  NLK 
(1H.6  kHz),  Jim  Creek,  Washington  to  Brisbane,  Australia.  The  variation  of 
the  geomagnetic  field  along  the  propagation  paths  are  shown  in  figure  17.  The 
individual  data  [xiints  on  the  plots  denote  locations  at  which  modal  constants 
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were  computed.  In  both  cases,  the  lower  boundary  is  assumed  to  be  sea  water 
with  a  conductivity  of  4.64  mho/m  and  a  relative  dielectric  constant  of  81. 

1 .  NWC  to  San  Diego 

a.  Nighttime 

The  propagation  constants  for  the  assumed  night  profile  for  the  NWC  path 

are  shown  in  figure  18.  Recall  that  the  x-variation  is  assumed  to  be  exp(- 

ikQSnx),  where  the  complex  propagation  factor  Sn  is  equal  to  sin0n,  and  0  is 

the  (complex)  eigenangle.  The  eigenangle  used  here  is  referred  to  the 

altitude  of  H  =  50  km  in  the  "transformed"  gui  'e  with  refractive  index  ez/a. 

From  Snell's  Law,  eH^a  sinO'  =  sinG„  where  O'  is  the  eigenangle'  discussed 

n  n  n 

further  in  this  study.  Figure  18a  shows  the  real  part  of  the  eigenangles 
(®r)»  while  figure  18b  shows  the  imaginary  part  (0^),  both  plotted  as  a 
function  of  distance  from  the  transmitter.  The  eigenangles  appear  to  fall  into 
two  groups.  One  group,  the  odd  numbered  modes,  is  cha racter ized  by  a  nearly 
constant  9r  with  perhaps  a  slight  increase  to  a  maximum  at  6  megameters  (the 
geomagnetic  equator).  The  for  these  odd  numbered  modes  diminishes  in  mag¬ 
nitude  as  the  magnetic  equator  is  approached.  The  even  numbered  modes  have 
just  the  opposite  behavior.  Note  that  the  mode  numbering  is  arbitrary;  it  is 
chosen  to  increase  monotonica 1 ly  with  increasing  0  at  the  transmitter. 

There  is  a  fundamental  physical  difference  between  the  even  and  odd 
numbered  modes  for  this  case.  This  is  seen  by  a  consideration  of  the  polari¬ 
zation  mixing  ratio,  defined  by  Pappert  (19f.8).  This  function  is  the  magni¬ 
tude  of  the  reciprocal  of  equation  (25)  and  is  approximately  equal  to  the 
ratio  of  the  maximum  magnitude  of  h  to  the  maximum  magnitude  of  e  within  the 

y  y 

guide.  Thus,  a  very  large  mixing  ratio  denotes  a  principally  TM  polarization 
whereas  a  very  small  ratio  denotes  principally  TE  pol a r izat ion.  The  mixing 
ratios  for  this  path  are  shown  in  figure  19.  The  even  numbered  modes  are  seen 
to  be  principally  TE,  dominantly  so  at  the  equator.  The  odd  numbered  modes 
are  principally  TM  polarized. 

The  propagation  characteristics  of  the  two  groups  of  modes  are  distinctly 
different.  Mode  constants  are  plotted  in  figures  20,  21  and  22.  For  the 

various  modes,  figure  20  shows  the  phase  velocities  and  figure  21  the  attenua¬ 
tion  rates.  Note  that  the  phase  velocities  and  attenuation  rates  can  also  be 
grouped  into  sets,  even  modes  and  odd  modes.  The  even  (TE)  modes  exhibit  a 
maximum  in  both  phase  velocity  and  attenuation  rate  at  the  equator  while  the 


odd  (TM)  modes  do  just  the  opposite.  The  increase  in  even  mode  phase  velocity 
indicates  a  decrease  in  the  effective  ionospheric  height.  Decreasing  the 
ionospheric  height  typically  results  in  an  increased  attenuation  rate,  and 
this  is  what  occurs.  Recall  from  Chapter  I  that  the  modal  excitation  factor 
is  conventionally  defined  for  a  vertical  dipole  source  exciting  a  principally 
vertically  polarized  mode.  Thus,  from  the  polarization  mixing  ratio,  the  odd 
modes  (TM)  can  be  expected  to  increase  in  excitation  efficiency  at  the  equator 
and  the  even  modes  can  be  expected  to  decrease.  This  is  indeed  the  case,  as 
is  indicated  in  figure  22. 

It  might  be  conjectured  from  the  excitation  factors  that  the  mode  sum  for 
this  case  would  exhibit  a  highly  moded  structure  over  the  entire  path 
length.  However,  examination  of  the  attenuation  rates  shows  that  the  highly 
excited  odd  modes  also  have  high  attenuation  rates  except  close  to  the  equa¬ 
tor.  Because  mode  1  has  a  relatively  small  excitation,  the  dominant  mode  is 
mode  3.  Mode  5  is  second  dominant  to  a  distance  of  about  10  megameters,  where 
mode  1  begins  to  increase  in  importance.  The  mode  sums  for  this  case  are 
presented  in  figure  23.  The  WKB  mode  sum  is  presented  along  with  the  mode  sum 
computed  using  the  formulation  of  Chapter  II.  The  two  sums  are  very  nearly 
identical,  indicating  essentially  no  mode  coupling  along  the  propagation 
path.  This  is  a  further  indication  of  the  correctness  of  the  computer  code 
used  in  this  study. 

b.  Sunrise/Sunset  Transition 

Numerical  modeling  of  transition  fading  requires  an  ionosphere  model  in 
the  neighborhood  of  the  terminator.  The  model  used  assumes  a  series  of 
discrete  steps  approximat ing  a  linear  variation  in  both  3  and  h' .  A  prelim¬ 
inary  investigation  of  transition  fading  was  made  by  considering  an  approxi¬ 
mate  path  with  constant  magnetic  conditions  and  a  variable  terminator  thick¬ 
ness.  As  a  consequence  of  this  approximation,  the  mode  conversion  coeffi¬ 
cients  for  a  given  terminator  thickness  are  independent  of  the  location  of  the 
terminator  along  the  path. 

Computed  amplitudes  and  phases  of  transition  fading  at  San  Diego  (15,000 
km  from  the  transmitter)  are  shown  in  figure  24  for  four  different  terminator 
dimensions.  Figure  24a  shows  the  sunrise  results  while  figure  24b  shows 
results  for  sunset.  Note  that  the  figures  have  the  distance  axes  reversed  so 
that  time  is  modeled  as  progressing  from  left  to  right.  For  both  sunrise  and 


sunset,  the  spacing  between  deep  fades  is  approximately  2000  km  for  all  values 
to  the  terminator  thickness.  The  absolute  position  of  each  fade,  however, 
does  depend  on  the  value  of  the  terminator  thickness.  There  is  approximately 
a  400  km  difference  in  the  fade  position  among  the  various  terminator 
thicknesses. 

The  fading  results  of  figure  24,  when  compared  with  the  mode  sums  of 
figure  23,  point  out  an  important  result.  The  mode  sum  shows  an  interference 
spacing  of  approximately  1000  km  over  most  of  the  path  length.  This  might 
lead  to  a  similar  expectation  for  the  transition  fade  spacing.  However,  the 
transition  fade  spacing  depends  primarily  on  the  difference  between  the  phase 
velocities  of  the  modes  that  are  most  efficiently  converted  by  the  terminator, 
and  not  on  all  the  modes  that  are  involved  in  the  pre-terminator  mode  sum. 
Examination  of  figure  23  reveals  the  existence  of  a  dominant  null  spacing  of 
2000  km  with  a  subdominant  mode  that  is  gradually  decaying.  Also,  while  the 
interference  nulls  for  the  mode  sum  are  about  5  dB  in  depth,  the  transition 
fades  are  almost  10  dB  in  depth  and  often  as  much  as  20  to  40  dB . 

The  change  of  phase  during  sunrise  shows  an  increase  over  the  first  half 
of  the  sunrise  transition  and  then  a  decrease  over  the  second  half.  This 
implies  that  the  phase  velocity  of  the  dominant  mode  was  modified  at  midtran¬ 
sition,  and  this  demonstrates  that  different  modes  were  dominant  in  the  early 
transition  time  and  in  the  late  transition  time.  The  sunset  phase  follows  a 
similar  behavior  except  for  the  first  fade  with  transition  thickness  of  750 
km,  where  a  phase  reversal  occurs. 

Signal  amplitudes  as  a  function  of  terminator  position  along  the  NWC  to 
San  Diego  propagation  path,  taking  into  account  variations  of  the  earth's 
magnetic  field,  have  been  confuted  for  several  distances  from  the  transmitter 
and  are  shown  in  figures  25  and  26.  The  transition  from  h'  =  74  km  to  h*  =  86 
km  (sunset)  is  shown  in  figure  25,  while  the  h'  =  86  km  to  h'  =  74  km  (sun¬ 
rise)  case  is  shown  in  figure  26.  The  terminator  width  was  taken  as  500  km  in 
both  cases.  The  zero  amplitude  for  each  successive  receiver  position  is 
offset  5  dB  from  that  for  the  previous  receiver  in  order  to  improve  visual 
clarity.  The  figures  support  completely  the  simplified  Crombie  model  of 
transition  fading.  In  figure  25,  the  fades  all  occur  when  the  receiver  is  a 
fixed  distance  from  the  terminator  whereas  in  figure  26  all  receiver  positions 
experience  a  fade  at  the  same  time  for  an  appropriate  position  of  the  termina— 
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tor.  Variations  from  these  patterns  occur,  of  course,  where  higher  order 
modes  are  important,  such  as  for  small  transmitter/terminator  distances  for 
sunset  and  for  large  distances  for  sunrise. 

The  modeling  results  for  the  NWC  to  San  Diego  propagation  path  (magneti¬ 
cally  northeasterly)  support  the  experimental  results  of  Chapter  III  in  all 
essential  details.  There  remains  now  to  apply  the  model  to  the  reverse  mag¬ 
netic  conditions,  that  is,  the  path  from  NLK  to  Brisbane. 

2.  NLK  to  Brisbane 

a.  Nighttime 

The  propagation  constants  for  the  assumed  night  profile  for  the  NLK  path 
are  presented  in  figure  27.  Figure  27a  shows  the  real  part  of  the  eigenangles 
(9r),  the  behavior  of  which  is  slightly  more  complicated  than  that  for  the  NWC 
to  San  Diego  path.  The  behavior  of  the  imaginary  part  of  the  eigenangles  (0^) 
in  figure  27b  is  considerably  more  complicated  than  that  for  the  northeasterly 
path.  There  does  not  appear  to  be  a  grouping  in  the  0r*  With  the  exception 
of  modes  9  and  10,  the  odd  numbered  modes  increase  to  a  maximum  in  0r  at  the 
equator  (7.5  megameters),  while  the  even  numbered  decrease  to  a  minimum.  Mode 
10  behaves  as  an  odd  mode  while  mode  9  behaves  as  an  even  mode.  Recall  that 
the  mode  numbering  is  established  by  ordering  0r  at  the  transmitter.  If  the 
0r  were  ordered  for  some  other  location,  such  as  the  end  point  of  the  path, 
the  numbering  of  modes  9  and  10  would  be  reversed.  The  important  feature  is 
the  grouping  of  the  modes,  not  the  numbering  within  the  groups.  Although  not 
immediately  obvious,  the  0^  also  fit  into  two  groups.  The  even  modes  have 
less  variation  than  the  odd  and  are  grouped  together  with  small  values,  except 
in  very  close  to  the  equator.  The  odd  modes  become  most  negative  on  each  side 
of  the  equator,  and  then  increase  rapidly  at  the  equator. 

The  variations  of  polarization  with  distance  along  the  path  are  shown  in 
figure  28.  The  behavior  of  this  parameter  is  also  markedly  different  from 
that  for  the  northeasterly  case.  The  modes  begin  at  the  transmitter  in  two 
groups.  Modes  1,  4,  6,  8  and  9  are  principally  TE  at  the  transmitter  while 
the  other  modes  are  TM.  All  of  the  TM  modes,  except  mode  3,  become  principal¬ 
ly  TE  just  prior  to  the  equator,  and  even  mode  3  increases  its  admixture  of  TE 
polarization.  The  TM  admixture  is  then  increased  for  all  the  modes  through 
the  equator  and,  except  for  mode  2,  the  modes  attain  a  maximum  of  TM  polariza- 
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tion  just  following  the  equator.  Then  the  modal  polarizations  return  essen¬ 
tially  to  the  conditions  at  the  transmitter.  The  behavior  of  the  modal 
polarization  exerts  a  strong  influence  on  the  mode  conversion  properties  of 
the  terminator. 

Because  the  daytime  ionosphere  is  only  slightly  anisotropic,  the  daytime 
modes  are  essentially  purely  TE  or  purely  TM.  Further,  a  purely  TE  mode  does 
not  convert  to  a  purely  TM  mode.  Thus,  the  conversion  of  the  essentially  TE 
nighttime  modes  at  about  7  megameters  from  the  transmitter  to  the  TM  daytime 
modes  can  be  expected  to  be  quite  different  from  the  conversion  of  the 
essentially  TM  nighttime  modes  at  about  8  megameters. 

The  modal  phase  velocities  and  attenuation  rates  are  shown  in  figures  29 
and  30.  The  phase  velocities  (figure  29)  exhibit  the  following  characteris¬ 
tic.  Not  only  do  the  phase  velocities  for  modes  9  and  10  cross  near  the 
transmitter  as  might  be  expected  from  the  0r  plots,  but  there  are  also  cros¬ 
sing  near  the  equator  for  modes  4,  5  and  6,  7.  A  mode  numbering  scheme  most 
commonly  favored  is  to  order  modes  according  to  phase  velocity,  smallest  to 
largest.  But  this  is  unfeasible  when  the  modal  phase  velocities  cross  at  the 
equator . 

Attenuation  rates  are  shown  in  figure  30.  Modes  1,  2,  and  4  exhibit 
similar  behavior  with  maximum  attenuation  rates  at  the  equator.  The  remaining 
modes  have  maxima  at  pre-and  post-equator  locations,  with  minima  at  the  equa¬ 
tor.  Mode  7  shows  a  slight  recovery  at  the  equator,  as  does  mode  10,  which 
has  an  attenuation  rate  larger  than  the  maximum  plotted  scale.  Clearly,  the 
attenuation  rate  has  no  apparent  correlation  to  the  polarization.  This  behav¬ 
ior  for  the  southwesterly  direction  is  in  complete  contrast  to  that  for  the 
northeasterly  direction. 

Mode  sums  for  our  nighttime  model  of  the  ionosphere  for  this  propagation 
path  are  shown  in  figure  31.  Both  the  WKB  and  the  mode  conversion  results  are 
presented.  In  contrast  to  the  behavior  for  the  northeasterly  path,  the  WKB 
and  mode  conversion  results  are  not  in  agreement  for  path  positions  near  and 
beyond  the  equator.  For  distances  shorter  than  about  6  megameters,  the  two 
mode  sums  are  in  complete  agreement  and  both  indicate  a  dominant  modal  inter¬ 
ference  spacing  of  approximtely  2  megameters.  This  spacing  is  not  apparent 
however,  for  distances  further  removed  from  the  transmitter. 
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b.  Sunrise/Sunset  Transition 


Signal  amplitudes  as  a  function  of  terminator  position  along  the  NLK  to 
Brisbane  path,  taking  into  account  variations  of  the  earth's  magnetic  field, 
have  been  completed  for  several  distances  from  the  transmitter  using  a 
terminator  of  500  km  width.  Let  us  first  consider  a  transition  from  h'  =  74 
to  h'  =86.  Because  the  path  direction  is  reversed  from  the  NWC  to  San  Diego 
case,  this  transition  corresponds  to  sunrise  rather  than  sunset.  Signal 
amplitudes  as  a  function  of  terminator  position  along  the  path  are  presented 
in  figure  32  for  several  receiver  locations.  The  remarkable  dissimilarity 
between  these  results  and  the  northwesterly  results  (figure  25)  is  immediately 
apparent.  A  number  of  points  should  be  noticed.  Receiver  distances  12  mega¬ 
meters  or  less  have  two  clearly  defined  fading  nulls,  whereas  receivers  much 
further  from  the  transmitter  apparently  exhibit  three  fading  nulls.  Further, 
the  double  null  spacings  are  different  for  each  receiver  location.  For 

terminator  positions  outside  the  equatorial  region  (±20°  latitude  as  shown), 
the  last  fade  spacing  for  the  triple  fade  cases  is  approximately  2  megameters 
while  the  fade  spacing  within  the  equatorial  region  is  nearly  3.5  megameters. 
All  of  these  observations  are  in  complete  agreement  with  the  experimental 
results  of  Lynn  (1969)  on  multiple  observations  of  the  equatorial  anomaly. 
Thus,  the  transequator ial  sunrise  transition  reported  by  Lynn  as  anomalous  are 
completely  explained  by  correctly  allowing  for  variation  of  the  geomagnetic 
field  along  the  path. 

The  sunset  results  ( h *  =  86  km  to  h'  =  74  km)  are  presented  in  figure 
33.  In  this  case,  there  are  two,  and  possibly  three,  fades  for  the  terminator 
on  the  transmitter  side  of  the  equatorial  region.  The  fade  spacing  is  approx¬ 
imately  2  megameters,  in  agreement  with  what  is  expected  for  midlatitude 
fading.  Further,  the  fades  are  simultaneous  at  all  receiver  sites,  also  in 
agreement  with  what  is  expected  for  midlatitudes.  When  the  terminator  is  in 
the  equatorial  region,  a  deep  fade  occurs  at  all  receiver  sites  nearly  simul¬ 
taneously.  This  fade  is  at  a  distance  of  approximately  3.5  megameters  of 
terminator  movement  from  the  previous  very  slight  fade.  These  results,  unfor¬ 
tunately,  do  not  correspond  with  those  of  Lynn  (1969),  who  reported  no  transi¬ 
tion  fading  at  all  for  the  sunset  case. 


V.  CONCLUSIONS 


This  study  reports  on  observations  of  sunrise/sunset  fading  on  a  long, 
transequator ial  northeasterly  propagation  path.  A  mode  conversion  model  is 
developed  and  numerical  modeling  results  are  presented  for  transition  fading 
on  long  transequator ial  paths,  both  northeasterly  and  southwesterly. 

From  these  results  the  following  general  conclusions  are  drawn: 

a.  Anomalous  transequator ial  transition  fading  is  not  observed  for  the 
northeasterly  propagation  path.  The  transition  fading  data  support  completely 
the  mode  conversion  model  of  transition  fading  developed  by  Crombie  (1964). 

b.  The  procedures  of  terminating  the  integration  limits  on  the  orthog¬ 
onal  i£at  ion  integrals  below  the  ionosphere  at  a  height  within  the  earth- 
ionosphere  waveguide  provides  a  considerable  simplification  for  numerical 
modeling  of  mode  conversion  effects.  Furthermore,  the  procedure  provides 
acceptably  accurate  results  for  such  modeling  studies. 

c.  The  transequator ial  sunrise  transitions  reported  in  the  literature  as 
anomalous  are  completely  explained  by  correctly  allowing  for  variations  of  the 
geomagnetic  field  along  the  path.  The  geomagnetic  field  influences  the  night¬ 
time  modal  parameters  to  such  an  extent  for  the  southwesterly  path  considered 
that  the  approximate  WKB  method  of  mode  summation  cannot  be  employed. 

d.  Careful  application  of  model  techniques  appears  to  be  capable  of 
describing  all  LF  and  VLF  propagation  in  the  earth- ionosphere  waveguide, 
provided  the  effect  of  the  earth's  magnetic  field  is  properly  included  at 
night.  However,  an  improved  model  of  the  ionosphere  at  altitudes  below  100  km 
is  much  to  be  desired. 
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Figure  23.  Mode  Sums  tor  NWC  to  San  Diego. 
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Figure  24b.  Computed  Transition  Fading  -  NWC  to  San  Diego  Day  to  Night  (h’  =  74  km  to  h'  =  Xf>  km). 
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Day  to  Night  (h'  =  74  km  to  h'  =  8(>  km). 
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POLARIZATION  MIXING  RATIO 
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Figure  30.  Attenuation  Rates  for  NLK  to  Brisbane 
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TRANSMITTER  TERMINATOR  DISTANCE  (MEGAMETERS) 


